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We simulate numerically Boussinesq convection in non-rotating spherical shells for a fluid 
with a unity Prandtl number and Rayleigh numbers up to 10®. In this geometry, curvature 
and radial variations of the gravitational acceleration yield asymmetric boundary layers. 
A systematic parameter study for various radius ratios (from r] = Vi/ro = 0.2 to 77 = 
0.95) and gravity profiles allows us to explore the dependence of the asymmetry on 
these parameters. We find that the average plume spacing is comparable between the 
spherical inner and outer bounding surfaces. An estimate of the average plume separation 
allows us to accurately predict the boundary layer asymmetry for the various spherical 
shell configurations explored here. The mean temperature and horizontal velocity profiles 
are in good agreement with classical Prandtl-Blasius laminar boundary layer profiles, 
provided the boundary layers are analysed in a dynamical frame that fluctuates with 
the local and instantaneous boundary layer thicknesses. The scaling properties of the 
Nusselt and Reynolds numbers are investigated by separating the bulk and boundary 
layer contributions to the thermal and viscous dissipation rates using numerical models 
with 77 = 0.6 and a gravity proportional to 1/r^. We show that our spherical models are 
consistent with the predictions of Grossmann & Lohse’s (2000) theory and that Nu[Ra) 
and Re{Ra) scalings are in good agreement with plane layer results. 
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1. Introduction 

Thermal convection is ubiquitous in geophysical and astrophysical fluid dynamics and 
rules, for example, turbulent flows in the interiors of planets and stars. The so-called 
Rayleigh-Benard (hereafter RB) convection is probably the simplest paradigm to study 
heat transport phenomena in these natural systems. In this configuration, convection is 
driven in a planar fluid layer cooled from above and heated from below (figure la). The 
fluid is confined between two rigid impenetrable walls maintained at constant tempera¬ 
tures. The key issue in RB convection is to understand the turbulent transport mecha¬ 
nisms of heat and momentum across the layer. In particular, how does the heat transport, 
characterised by the Nusselt number Nu, and the flow amplitude, characterised by the 
Reynolds number Re, depend on the various control parameters of the system, namely 
the Rayleigh number Ra, the Prandtl number Pr and the cartesian aspect ratio T? In 
general, T = W/H quantifies the fluid layer width W over its height H in classical planar 
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Figure 1. Schematic showing different Rayleigh-Benard convection setups, (a) An infinitely 
extended fluid layer of height H heated from below and cooled from above, (b) The typical RB 
convection setup in a cylindrical cell of aspect ratio F = W/H. (c) Convection in a spherical 
shell with a radius ratio r] = Ti/ro with a radially inward gravity. In the three panels, the red 
and the blue surfaces correspond to the hot and cold boundaries held at constant temperatures. 


or cylindrical RB cells. In spherical shells, we rather employ the ratio of the inner to the 
outer radius rj = Ti/ro to characterise the geometry of the fluid layer. 

Laboratory experiments of RB convection are classically performed in rectangular or in 
cylindrical tanks with planar upper and lower bounding surfaces where the temperature 
contrast is imposed (see figure 16). In such a system, the global dynamics are strongly 
influenced by the flow properties in the thermal and kinematic boundary layers that form 
in the vicinity of the walls. The characterisation of the structure of these boundary layers 
is crucial for a better understanding of the transport processes. The marginal stability 
theory by Malkus (1954) is the earliest boundary layer model and relies on the assumption 
that the thermal boundary layers adapt their thicknesses to maintain a critical boundary 
layer Rayleigh number, which implies Nu ^ Ra}/^. Assuming that the boundary layers 
are sheared, Shraiman & Siggia (1990) later derived a theoretical model that yields 
scalings of the form Nu ~ Ra}!’^ and Re ^ (see also Siggia 1994). 

These asymptotic laws were generally consistent with most of the experimental results 
obtained in the 1990s up to Ra < 10^^. Within the typical experimental resolution of 
one percent, simple power laws of the form Nu ^ RaPPr^ were found to provide an 
adequate representation with a exponents ranging from 0.28 to 0.31, in relatively good 
agreement with the Shraiman & Siggia model (e.g. Castaing et al. 1989; Chavanne et al. 
1997; Niemela et al. 2000). However, later high-precision experiments by Xu et al. (2000) 
revealed that the dependence of Nu upon Ra cannot be accurately described by such 
simple power laws. In particular, the local slope of the function Nu{Ra) has been found 
to increase slowly with Ra. The effective exponent aes of Nu ^ Ra°^'‘^^ roughly ranges 
from values close to 0.28 near Ra ~ 10^ — 10® to 0.33 when Ra ~ 10^^ — 10^^ (e.g. 
Funfschilling et al. 2005; Cheng et al. 2015). 

Grossmann & Lohse (2000, 2004) derived a competing theory capable of capturing 
this complex dynamics (hereafter GL). This scaling theory is built on the assumption 
of laminar boundary layers of Prandtl-Blasius (PB) type (Prandtl 1905; Blasius 1908). 
According to the GL theory, the flows are classified in four different regimes in the Ra—Pr 
phase space according to the relative contribution of the bulk and boundary layer viscous 
and thermal dissipation rates. The theory predicts non-power-law behaviours for Nu and 
Ra in good agreement with the dependence Nu = f{Ra,Pr,T) and Re = f{Ra,Pr,r) 
observed in recent experiments and numerical simulations of RB convection in planar 
or cylindrical geometry (see for recent reviews Ahlers et al. 2009; Chilla & Schumacher 
2012 ). 

Benefiting from the interplay between experiments and direct numerical simulations 
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(DNS), turbulent RB convection in planar and cylindrical cells has received a lot of 
interest in the past two decades. However, the actual geometry of several fundamental 
astrophysical and geophysical flows is essentially three-dimensional within concentric 
spherical upper and lower bounding surfaces under the influence of a radial buoyancy 
force that strongly depends on radius. The direct applicability of the results derived in 
the planar geometry to spherical shell convection is thus questionable. 

As shown in figure 1(c), convection in spherical shells mainly differs from the tradi¬ 
tional plane layer configuration because of the introduction of curvature and the absence 
of side walls. These specific features of thermal convection in spherical shells yield sig¬ 
nificant dynamical differences with plane layers. For instance, the heat flux conservation 
through spherical surfaces implies that the temperature gradient is larger at the lower 
boundary than at the upper one to compensate for the smaller area of the bottom sur¬ 
face. This yields a much larger temperature drop at the inner boundary than at the outer 
one. In addition, this pronounced asymmetry in the temperature profile is accompanied 
by a difference between the thicknesses of the inner and the outer thermal boundary 
layers. Following Malkus’s marginal stability arguments, Jarvis (1993) and Vangelov & 
Jarvis (1994) hypothesised that the thermal boundary layers in curvilinear geometries 
adjust their thickness to maintain the same critical boundary layer Rayleigh number 
at both boundaries. This criterion is however in poor agreement with the results from 
numerical models (e.g. Deschamps et al. 2010). The exact dependence of the boundary 
layer asymmetry on the radius ratio and the gravity distribution thus remains an open 
question in thermal convection in spherical shells (Bercovici et al. 1989; Jarvis et al. 1995; 
Sotin & Labrosse 1999; Shahnas et al. 2008; O’Farrell et al. 2013). This open issue sheds 
some light on the possible dynamical influence of asymmetries between the hot and cold 
surfaces that originate due to both the boundary curvature and the radial dependence 
of buoyancy in spherical shells. 

Ground-based laboratory experiments involving spherical geometry and a radial buoy¬ 
ancy forcing are limited by the fact that gravity is vertically downwards instead of radially 
inwards (Scanlan et al. 1970; Feldman & Colonius 2013). A possible way to circumvent 
this limitation is to conduct experiments under microgravity to suppress the vertically 
downward buoyancy force. Such an experiment was realised by Hart et al. (1986) who 
designed a hemispherical shell that flew on board of the space shuttle Challenger in May 
1985. The radial buoyancy force was modelled by imposing an electric field across the 
shell. The temperature dependence of the fluid’s dielectric properties then produced an 
effective radial gravity that decreases with the fifth power of the radius (i.e. g ~ l/''”^)- 
More recently, a similar experiment named “GeoFlow” was run on the International 
Space Station, where much longer flight times are possible (Futterer et al. 2010; Fut- 
terer et al. 2013). This later experiment was designed to mimic the physical conditions 
in the Earth mantle. It was therefore mainly dedicated to the observation of plume-like 
structures in a high Prandtl number regime {Pr > 40) for Ra < 10®. Unfortunately, this 
limitation to relatively small Rayleigh numbers makes the GeoFlow experiment quite 
restricted regarding asymptotic scaling behaviours in spherical shells. 

To compensate the lack of laboratory experiments, three dimensional numerical models 
of convection in spherical shells have been developed since the 1980s (e.g. Zebib et al. 
1980; Bercovici et al. 1989, 1992; Jarvis et al. 1995; Tilgner 1996; Tilgner & Busse 1997; 
King et al. 2010; Choblet 2012). The vast majority of the numerical models of non¬ 
rotating convection in spherical shells has been developed with Earth’s mantle in mind. 
These models therefore assume an infinite Prandtl number and most of them further 
include a strong dependence of viscosity on temperature to mimic the complex rheology of 
the mantle. Several recent studies of isoviscous convection with infinite Prandtl number in 
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spherical shells have nevertheless been dedicated to the analysis of the scaling properties 
of the Nusselt number. For instance, Deschamps et al. (2010) measured convective heat 
transfer in various radius ratios ranging from 77 = 0.3 to rj = 0.8 and reported Nu ~ 
RaO.273 2 q4 < while Wolstencroft et al. (2009) computed numerical models 

with Earth’s mantle geometry (77 = 0.55) up to Ra = 10® and found Nu ^ 

These studies also checked the possible influence of internal heating and reported quite 
similar scalings. 

Most of the numerical models of convection in spherical shells have thus focused on the 
very specific dynamical regime of the infinite Prandtl number. The most recent attempt 
to derive the scaling properties of Nu and Re in non-rotating spherical shells with finite 
Prandtl numbers is the study of Tilgner (1996). He studied convection in self-graviting 
spherical shells (i.e. g ^ r) with 77 = 0.4 spanning the range 0.06 < Pr < 10 and 
4 X 10® < Ra < 8 X 10®. This study was thus restricted to low Rayleigh numbers, 
relatively close to the onset of convection, which prevents the derivation of asymptotic 
scalings for Nu{Ra, Pr) and Re{Ra, Pr) in spherical shells. 

The objectives of the present work are twofold: (i) to study the scaling properties of 
Nu and Re in spherical shells with finite Prandtl number; (ii) to better characterise the 
inherent asymmetric boundary layers in thermal convection in spherical shells. We there¬ 
fore conduct two systematic parameter studies of turbulent RB convection in spherical 
shells with Pr = 1 by means of three dimensional DNS. In the first set of models, we vary 
both the radius ratio (from 77 = 0.2 to 77 = 0.95) and the radial gravity profile (consider¬ 
ing g(r) G [T’/ro, 1 , (rg/r)'^, (t'o/t')®]) in a moderate parameter regime (i.e. 5 < Nu < 15 
for the majority of the cases) to study the influence of these properties on the boundary 
layer asymmetry. We then consider a second set of models with 77 = 0.6 and g ^ 1/r^ 
up to Ra = 10®. These DNS are used to check the applicability of the GL theory to 
thermal convection in spherical shells. We therefore numerically test the different basic 
prerequisites of the GL theory: we first analyse the nature of the boundary layers before 
deriving the individual scaling properties for the different contributions to the viscous 
and thermal dissipation rates. 

The paper is organised as follows. In § 2, we present the governing equations and the 
numerical models. We then focus on the asymmetry of the thermal boundary layers in 
§ 3. In § 4, we analyse the nature of the boundary layers and show that the boundary 
layer profiles are in agreement with the Prandtl-Blasius theory (Prandtl 1905; Blasius 
1908). In § 5, we investigate the scaling properties of the viscous and thermal dissipation 
rates before calculating the Nu{Ra) and Re{Ra) scalings in § 6. We conclude with a 
summary of our findings in § 7. 


2. Model formulation 

2.1. Governing hydrodynamical equations 

We consider RB convection of a Boussinesq fluid contained in a spherical shell of outer 
radius ro and inner radius r.^. The boundaries are impermeable, no slip and at constant 
temperatures T\,ot and Ttop- We adopt a dimensionless formulation using the shell gap 
d = ro — ri as the reference lengthscale and the viscous dissipation time d^/v as the 
reference timescale. Temperature is given in units of AT = Ttop — T},ot, the imposed 
temperature contrast over the shell. Velocity and pressure are expressed in units oivjd 
and poV^/d^, respectively. Gravity is non-dimensionalised using its reference value at the 
outer boundary go- The dimensionless equations for the velocity u, the pressure p and 
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the temperature T are given by 


V • M = 0, 


„ Ra ^ 

— +u-Vu = -Vp+ —gTer 
at Fr 

dT 1 

-Fu.VT=-AT, 


Am, 


( 2 . 1 ) 

( 2 . 2 ) 


(2.3) 


where is the unit vector in the radial direction and g is the gravity. Several gravity 
profiles have been classically considered to model convection in spherical shells. For in¬ 
stance, self-graviting spherical shells with a constant density correspond to g r (e.g 
Tilgner 1996), while RB convection models with infinite Prandtl number usually assume 
a constant gravity in the perspective of modelling Earth’s mantle (e.g. Bercovici et al. 
1989). The assumption of a centrally-condensed mass has also been frequently assumed 
when modelling rotating convection (e.g. Gilman & Glatzmaier 1981; Jones et al. 2011) 
and yields g ^ 1/r^. Finally, the artificial central force field of the microgravity experi¬ 
ments takes effectively the form oi g ^ 1/r® (Hart et al. 1986; Feudel et al. 2011; Futterer 
et al. 2013). To explore the possible impact of these various radial distribution of buoy¬ 
ancy on RB convection in spherical shells, we consider different models with the four 
following gravity profiles: g € [r/ro, 1, (ro/r)^, (ro/r)®]. Particular attention will be paid 
to the cases with g = (ro/r)^, which is the only radial function compatible with an exact 
analysis of the dissipation rates (see below, § 2.3). 

The dimensionless set of equations (2. 1-2.3) is governed by the Rayleigh number Ra, 
the Prandtl number Pr and the radius ratio of the spherical shell rj defined by 


Ra = 


ag. 


,ATd^ 


Pr = -, 

K 


rj = 


(2.4) 


where ir and k are the viscous and thermal diffusivities and a is the thermal expansivity. 


2.2. Diagnostic parameters 

To quantify the impact of the different control parameters on the transport of heat and 
momentum, we analyse several diagnostic properties. We adopt the following notations 
regarding different averaging procedures. Overbars correspond to a time average 


/ 


1 

T 



/dt. 


where r is the time averaging interval. Spatial averaging over the whole volume of the 
spherical shell are denoted by triangular brackets (•••), while {■■■)s correspond to an 
average over a spherical surface: 


if) = ^ = /(r,6»,(())sin6»d6»d(^. 


where V is the volume of the spherical shell, r is the radius, 6 the colatitude and (j) the 
longitude. 

The convective heat transport is characterised by the Nusselt number Nu, the ratio of 
the total heat flux to the heat carried by conduction. In spherical shells with isothermal 
boundaries, the conductive temperature profile Tc is the solution of 


_d 

dr 



= 0 , 


Tc{r^) = 1 , 


Tc{ro) = 0 , 
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which yields 


Tc{r) = 


V 


V 




1 - ry 


(2.5) 


For the sake of clarity, we adopt in the following the notation d for the time-averaged 
and horizontally-averaged radial temperature profile: 


d(r) = (T), 


The Nusselt number then reads 


Nu = 


1 dT, 

Pr dr 


=_Ar = n) = —^{T = T,). 
dr Tj dr 


( 2 . 6 ) 


The typical rms flow velocity is given by the Reynolds number 


Re = ^/{vA) = yl{u^ + ul+ ul), (2.7) 

while the radial profile for the time and horizontally-averaged horizontal velocity is de¬ 
fined by 

Rehir) = + ( 2 - 8 ) 


2.3. Exact dissipation relationships in spherical shells 


The mean buoyancy power averaged over the whole volume of a spherical shell is ex¬ 
pressed by 


P = 


Ra 

Pr 


{g UrT) 


Att Ra 
V Pr 


gr"^ (urT),dr, 


Using the Nusselt number definition (2.6) and the conductive temperature profile (2.5) 
then leads to 


47 r Ra 

V Pr^ 


gr 


— dr — Nu 
dr 


V 

{1-gy 



The first term in the parentheses becomes identical to the imposed temperature drop 
AT for a gravity g ^ 1/r^: 



di? 

dr 


dr = rl [•d(ro) - d(ri)] 


—r 


2 

O J 


and thus yields an analytical relation between P and Nu. For any other gravity model, 
we have to consider the actual spherically-symmetric radial temperature profile 'd(r). 
Christensen & Aubert (2006) solve this problem by approximating d(r) by the diffusive 
solution (2.5) and obtain an approximate relation between P and -^{Nu — 1). This 
motivates our particular focus on the g = (j’o/r)'^ cases which allows us to conduct an 
exact analysis of the dissipation rates and therefore check the applicability of the GL 
theory to convection in spherical shells. 

Noting that /r ° 5 'dr = for g = {rolr)"^, one finally obtains the exact 

relation for the viscous dissipation rate eu- 

= = (2.9) 

The thermal dissipation rate can be obtained by multiplying the temperature equation 
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(2.3) by T and integrate it over the whole volume of the spherical shell. This yields 

eT = Wm = TT^^^Nu. ( 2 . 10 ) 

These two exact relations (2.9-2.10) can be used to validate the spatial resolutions of the 
numerical models with g = (ro/r)^. To do so, we introduce Xeu Xeti f^e ratios of 
the two sides of Eqs (2.9-2.10): 


{1 rj rj^) Pr^ 

^ 3i?a {Nu-l) ’ 
_ (1-^77-^772) 

377 iVu 


( 2 . 11 ) 


2.4. Setting up a parameter study 

2.4.1. Numerical technique 

The numerical simulations have been carried out with the magnetohydrodynamics 
code Magic (Wicht 2002). MagIC has been validated via several benchmark tests for 
convection and dynamo action (Christensen et al. 2001; Jones et al. 2011). To solve the 
system of equations (2. 1-2.3), the solenoidal velocity field is decomposed into a poloidal 
and a toroidal contribution 


M = V X (V X Wer) -I- V X Zer^ 

where W and Z are the poloidal and toroidal potentials. W, Z, p and T are then expanded 
in spherical harmonic functions up to degree £max in the angular variables 9 and (j) and in 
Chebyshev polynomials up to degree Nr in the radial direction. The combined equations 
governing W and p are obtained by taking the radial component and the horizontal part 
of the divergence of (2.2). The equation for Z is obtained by taking the radial component 
of the curl of (2.2). The equations are time-stepped by advancing the nonlinear terms 
using an explicit second-order Adams-Bashforth scheme, while the linear terms are time- 
advanced using an implicit Crank-Nicolson algorithm. At each time step, all the nonlinear 
products are calculated in the physical space (r, 9, (j)) and transformed back into the 
spectral space (r, £, m). For more detailed descriptions of the numerical method and 
the associated spectral transforms, the reader is referred to (Gilman & Glatzmaier 1981; 
Tilgner & Busse 1997; Christensen & Wicht 2007). 

2.4.2. Parameter choices 

One of the main focuses of this study is to investigate the global scaling properties of 
Pr = 1 RB convection in spherical shell geometries. This is achieved via measurements 
of the Nusselt and Reynolds numbers. In particular, we aim to test the applicability of 
the GL theory to spherical shells. As demonstrated before, only the particular choice of 
a gravity profile of the form g ^ l/r^ allows the exactness of the relation (2.9). Our 
main set of simulations is thus built assuming g = (toIv')^- The radius ratio is kept 
to 77 = 0.6 and the Prandtl number to Pr = 1 to allow future comparisons with the 
rotating convection models by Gastine & Wicht (2012) and Gastine et al. (2013) who 
adopted the same configuration. We consider 35 numerical cases spanning the range 
1.5 X 10^ < Ra < 10®. Table 1 summarises the main diagnostic quantities for this 
dataset of numerical simulations and shows that our models basically lie within the 
ranges 1 < i?e < 7000 and 1 < Nu < 70. 

Another important issue in convection in spherical shells concerns the determination 
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Ra 

Nu 

Re 

Xr/\°T 

Kj/>^u 

e ? 7 (%) 



X^U 

Nr X £rnax 

1.5 X 10 ® 

1.33 

4.4 

_ 

_ 

_ 

_ 

1.000 

1.000 

49 X 85 

2 X 10 ® 

1.59 

6.7 

- 

- 

- 

- 

1.000 

1.000 

49 X 85 

3 X 10 ® 

1.80 

9.6 

- 

- 

- 

- 

1.000 

1.000 

49 X 85 

5 X 10 ® 

2.13 

14.4 

- 

- 

- 

- 

1.000 

1.000 

49 X 85 

7 X 10 ® 

2.20 

17.5 

0 . 186 / 0.251 

0 . 076 / 0.104 

0.11 

0.57 

1.000 

1.000 

49 X 85 

9 X 10 ® 

2.43 

21.7 

0 . 168 / 0.223 

0 . 070 / 0.094 

0.14 

0.60 

1.000 

1.000 

49 X 85 

1 X 10 '‘ 

2.51 

23.3 

0 . 162 / 0.217 

0 . 069 / 0.092 

0.15 

0.61 

1.000 

1.000 

49 X 85 

1.5 X 10 ^ 

2.81 

29.8 

0 . 143 / 0.196 

0 . 062 / 0.086 

0.17 

0.62 

1.000 

1.000 

49 X 85 

2 X 10 '‘ 

3.05 

35.0 

0 . 130 / 0.185 

0 . 059 / 0.082 

0.15 

0.64 

1.000 

1.000 

49 X 85 

3 X 10 '‘ 

3.40 

44.0 

0 . 116 / 0.167 

0 . 054 / 0.077 

0.17 

0.64 

1.000 

1.000 

49 X 85 

5 X lO '^ 

3.89 

57.5 

0 . 102 / 0.147 

0 . 049 / 0.069 

0.18 

0.67 

1.000 

1.000 

49 X 85 

7 X lO '^ 

4.27 

68.5 

0 . 093 / 0.133 

0 . 046 / 0.062 

0.18 

0.69 

1.000 

1.000 

49 X 85 

1 X 10 ® 

4.71 

82.3 

0 . 085 / 0.120 

0 . 043 / 0.058 

0.18 

0.70 

1.000 

1.000 

49 X 85 

1.5 X 10 ® 

5.28 

101.2 

0 . 076 / 0.107 

0 . 039 / 0.053 

0.19 

0.71 

1.000 

1.000 

49 X 85 

2 X 10 ® 

5.72 

117.0 

0 . 070 / 0.099 

0 . 037 / 0.050 

0.20 

0.74 

1.000 

1.000 

49 X 85 

3 X 10 ® 

6.40 

143.3 

0 . 062 / 0.088 

0 . 033 / 0.046 

0.22 

0.74 

1.000 

1.000 

61 X 106 

5 X 10 ® 

7.37 

185.1 

0 . 054 / 0.077 

0 . 030 / 0.042 

0.24 

0.77 

1.000 

1.000 

61 X 106 

7 X 10 ® 

8.10 

218.6 

0 . 049 / 0.070 

0 . 028 / 0.040 

0.22 

0.77 

1.000 

1.000 

61 X 106 

1 X 10 ® 

8.90 

259.2 

0 . 045 / 0.064 

0 . 026 / 0.037 

0.24 

0.79 

1.000 

1.000 

81 X 170 

1.5 X 10 ® 

9.97 

315.1 

0 . 040 / 0.057 

0 . 024 / 0.034 

0.23 

0.81 

1.000 

1.000 

81 X 170 

2 X 10 ® 

10.79 

362.8 

0 . 037 / 0.053 

0 . 023 / 0.032 

0.26 

0.81 

1.000 

1.000 

81 X 170 

3 X 10 ® 

12.11 

443.5 

0 . 033 / 0.048 

0 . 020 / 0.030 

0.24 

0.82 

0.999 

1.003 

81 X 170 

5 X 10 ® 

13.97 

565.6 

0 . 029 / 0.041 

0 . 018 / 0.027 

0.25 

0.83 

1.000 

1.001 

97 X 266 

7 X 10 ® 

15.39 

666.4 

0 . 026 / 0.037 

0 . 017 / 0.025 

0.27 

0.83 

1.000 

1.005 

97 X 266 

1 X 10 ^ 

17.07 

790.4 

0 . 023 / 0.034 

0 . 016 / 0.024 

0.25 

0.84 

1.000 

1.005 

97 X 266 

1.5 X 10 '^ 

19.17 

960.1 

0 . 021 / 0.030 

0 . 015 / 0.021 

0.28 

0.84 

1.000 

1.009 

97 X 341 

2 X 10 '^ 

20.87 

1099.7 

0 . 019 / 0.028 

0 . 013 / 0.020 

0.27 

0.85 

1.000 

1.005 

121 X 426 

3 X lO ’’ 

23.50 

1335.5 

0 . 017 / 0.025 

0 . 012 / 0.018 

0.28 

0.85 

1.000 

1.012 

121 X 426 

5 X lO ’’ 

27.35 

1690.2 

0 . 014 / 0.021 

0 . 011 / 0.016 

0.27 

0.86 

1.000 

1.010 

161 X 512 

7 X 10 ^ 

30.21 

1999.1 

0 . 013 / 0.019 

0 . 010 / 0.015 

0.28 

0.86 

1.000 

1.005 

201 X 576 

1 X 10 ® 

33.54 

2329.9 

0 . 012 / 0.017 

0 . 009 / 0.013 

0.29 

0.87 

1.000 

1.011 

201 X 682 

2 X 10 ® 

41.49 

3239.3 

0 . 010 / 0.014 

0 . 008 / 0.012 

0.29 

0.88 

1.000 

1.006 

321 X 682 

3 X 10 ® 

47.22 

3882.5 

0 . 008 / 0.012 

0 . 007 / 0.010 

0.29 

0.87 

1.001 

1.015 

321 X 768 

5 X Id® 

56.67 

5040.0 

0 . 007 / 0.009 

0 . 005 / 0.008 

0.31 

0.85 

1.001 

1.050 

321 X 682 * 

5 X 10 ® 

55.07 

4944.1 

0 . 007 / 0.011 

0 . 006 / 0.009 

0.29 

0.88 

0.999 

1.003 

401 X 853 * 

1 X 10® 

73.50 

7039.4 

0 . 005 / 0.007 

0 . 004 / 0.006 

0.32 

0.79 

1.002 

1.148 

401 X 682 * 

1 X 10 ® 

68.48 

6802.5 

0 . 006 / 0.009 

0 . 005 / 0.007 

0.30 

0.89 

0.996 

1.006 

513 X 1066 * 


Table 1. Summary table of Pr = 1 numerical simulations with rj = 0.6 and g — (ro/r)^- 
The boundary layer thicknesses and the viscous and thermal dissipations are only given for 
the cases with Ra > 7 x 10^ when boundary layers can be clearly identihed. The italic lines 
indicate simulations with smaller truncations to highlight the possible resolution problems. The 
cases with Ra = 5 x 10® and Ra = 10® (highlighted with a star in the last column) have been 
computed assuming a two-fold symmetry, i.e. effectively simulating only half of the spherical 
shell in longitude. 


of the average bulk temperature and the possible boundary layer asymmetry between the 
inner and the outer boundaries (e.g. Jarvis 1993; Tilgner 1996). To better understand the 
influence of curvature and the radial distribution of buoyancy, we thus compute a second 
set of numerical models. This additional dataset consists of 113 additional simulations 
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with various radius ratios and gravity profiles, spanning the range 0.2 < rj < 0.95 with 
g S [r/r-Q, 1, (ro/r)^, (ro/r)®]. To limit the numerical cost of this second dataset, these 
cases have been run at moderate Rayleigh number and typically span the range 5 < 
Nu < 15 for the majority of the cases. Table 2, given in the Appendix, summarises the 
main diagnostic quantities for this second dataset of numerical simulations. 

2.4.3. Resolution checks 

Attention must be paid to the numerical resolutions of the DNS of RB convection (e.g. 
Shishkina et al. 2010). Especially, underresolving the fine structure of the turbulent flow 
leads to an overestimate of the Nusselt number, which then falsifies the possible scaling 
analysis (Amati et al. 2005). One of the most reliable ways to validate the truncations 
employed in our numerical models consists of comparing the obtained viscous and thermal 
dissipation rates with the average Nusselt number (Stevens et al. 2010; Lakkaraju et al. 
2012; King et al. 2012). The ratios x^u Xer; defined in (2.11), are found to be very 
close to unity for all the cases of Table 1, which supports the adequacy of the employed 
numerical resolutions. To further highlight the possible impact of inadequate spatial 
resolutions, two underresolved numerical models for the two highest Rayleigh numbers 
have also been included in Table 1 (lines in italics). Because of the insufficient number of 
grid points in the boundary layers, the viscous dissipation rates are significantly higher 
than expected in the statistically stationary state. This leads to overestimated Nusselt 
numbers by similar percent differences (2 — 10%). 

Table 1 shows that the typical resolutions span the range from {Nr = 49, £max = 85) 
to {Nr = 513, £niax = 1066). The two highest Rayleigh numbers have been computed 
assuming a two-fold azimuthal symmetry to ease the numerical computations. A com¬ 
parison of test runs with or without the two-fold azimuthal symmetry at lower Rayleigh 
numbers (5 x 10^ < Ra < 3 x 10®) showed no significant statistical differences. This en¬ 
forced symmetry is thus not considered to be influential. The total computational time 
for the two datasets of numerical models represents roughly 5 million Intel Ivy Bridge 
CPU hours. 


3. Asymmetric boundary layers in spherical shells 

3.1. Definitions 

Several different approaches are traditionally considered to define the thermal boundary 
layer thickness At. They either rely on the horizontally-averaged mean radial temperature 
profile 'd{r) or on the temperature fluctuation a defined as 



Among the possible estimates based on 'd(r), the slope method (e.g. Verzicco & Camussi 
1999; Breuer et al. 2004; Liu & Ecke 2011) defines At as the depth where the linear fit 
to i?(r) near the boundaries intersects the linear fit to the temperature profile at mid¬ 
depth. Alternatively, a exhibits sharp local maxima close to the walls. The radial distance 
separating those peaks from the corresponding nearest boundary can be used to define 
the thermal boundary layer thicknesses (e.g. Tilgner 1996; King et al. 2013). Figure 2(a) 
shows that both definitions of At actually yield nearly indistinguishable boundary layer 
thicknesses. We therefore adopt the slope method to define the thermal boundary layers. 

There are also several ways to define the viscous boundary layers. Figure 2(&) shows 
the vertical profile of the root-mean-square horizontal velocity Rch. This profile exhibits 
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Figure 2. (a) Radial profiles of the time and horizontally-averaged mean temperature i9(r) 
(solid black line) and the temperature variance a (dotted black line). The thermal boundary 
layers Ay and Ay are highlighted by the gray shaded area. They are defined as the depths where 
the linear fit to ■d(r) near the top (bottom) crosses the linear fit to the temperature profile 
at mid-depth (dashed black lines), (b) Radial profiles of the time and horizontally-averaged 
horizontal velocity Rch (r) . The viscous boundary layers are either defined by the local maxima 
of Reh (black dotted lines, ^u,m) ot by the intersection of the linear fit to Ren near the 

inner (outer) boundary with the tangent to the local maxima (dashed black lines). This second 
definition is highlighted by a gray shaded area (A)/, A^). These profiles have been obtained from 
a numerical model with Ra — 10^, rj — 0.6 and g — {rolr)^ ■ 


strong increases close to the boundaries that are accompanied by well-defined peaks. Fol¬ 
lowing Tilgner (1996) and Kerr & Herring (2000), the first way to define the kinematic 
boundary layer is thus to measure the distance between the walls and these local maxima. 
This commonly-used definition gives (.^u m) inner (outer) spherical bound¬ 

ary. Another possible method to estimate the viscous boundary layer follows a similar 
strategy as the slope method that we adopted for the thermal boundary layers (Breuer 
et al. 2004). A^ (A^) is defined as the distance from the inner (outer) wall where the 
linear fit to Reh near the inner (outer) boundary intersects the horizontal line passing 
through the maximum horizontal velocity. 

Figure 2{b) reveals that these two definitions lead to very distinct viscous boundary 
layer thicknesses. In particular, the definition based on the position of the local maxima 
of Reh yields much thicker boundary layers than the tangent intersection method, i.e. 
Kfm > ^u°- The discrepancies of these two definitions are further discussed in § 4. 

3.2. Asymmetric thermal boundary layers and mean bulk temperature 

Figure 2 also reveals a pronounced asymmetry in the mean temperature proHles with a 
much larger temperature drop at the inner boundary than at the outer boundary. As 
a consequence, the mean temperature of the spherical shell Tm = y dV is much 
below AT 12. Determining how the mean temperature depends on the radius ratio rj 
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Figure 3. (a) ^{r) for different radius ratios g with a gravity g = (vo/r)^. These models have 
approximately the same Nusselt number 12 < Nu < 14. (6) il(r) for different gravity proHles 
with a fixed radius ratio g = 0.6. These models have approximately the same Nusselt number 
10 < iVu < 11. 


has been an ongoing open question in mantle convection studies with infinite Prandtl 
number (e.g. Bercovici et al. 1989; Jarvis 1993; Vangelov & Jarvis 1994; Jarvis et al. 
1995; Sotin & Labrosse 1999; Shahnas et al. 2008; Deschamps et al. 2010; O’Farrell et al. 
2013). To analyse this issue in numerical models with Pr = 1, we have performed a 
systematic parameter study varying both the radius ratio of the spherical shell g and 
the gravity profile g{r) (see Table 2). Figure 3 shows some selected radial profiles of 
the mean temperature d for various radius ratios (panel a) and gravity profiles (panel 
b) for cases with similar Nu. For small values of g, the large difference between the 
inner and the outer surfaces lead to a strong asymmetry in the temperature distribution: 
nearly 90% of the total temperature drop occurs at the inner boundary when g = 0.2. In 
thinner spherical shells, the mean temperature gradually approaches a more symmetric 
distribution to finally reach Tm = 0.5 when 77 —>■ 1 (no curvature). Figure 3(6) also reveals 
that a change in the gravity profile has a direct impact on the mean temperature profile. 
This shows that both the shell geometry and the radial distribution of buoyancy affect 
the temperature of the fluid bulk in RB convection in spherical shells. 

To analytically access the asymmetries in thickness and temperature drop observed in 
figure 3, we first assume that the heat is purely transported by conduction in the thin 
thermal boundary layers. The heat flux conservation through spherical surfaces (2.6) 
then yields 


AT° 


.AT* 


= d 


A6 


(3.2) 


where the thermal boundary layers are assumed to correspond to a linear conduction 
profile with a temperature drop AT* (AT°) over a thickness (Ay). As shown in 
Figs. 2-3, the fluid bulk is isothermal and forms the majority of the fluid by volume. We 
can thus further assume that the temperature drops occur only in the thin boundary 
layers, which leads to 


AT° -b AT* = 1. 


(3.3) 


Equations (3.2) and (3.3) are nevertheless not sufficient to determine the three unknowns 
AT*, AT°, Ay/Ay and an additional physical assumption is required. 

A hypothesis frequently used in mantle convection models with infinite Prandtl number 
in spherical geometry (Jarvis 1993; Vangelov & Jarvis 1994) is to further assume that 
both thermal boundary layers are marginally stable such that the local boundary layer 
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Figure 4. (a) Ratio of boundary layer Rayleigh numbers (3.4) for various radius ratios and 
gravity profiles. The horizontal dashed line corresponds to the hypothetical identity Ra\ = Ral- 
(b) Temperature drop at the outer boundary layer. The lines correspond to the theoretical 
prediction given in (3.5). 


Rayleigh numbers Ra\ and Ra^ are equal: 


Ra\ = Ral 


agiAT'^Xj 


agoAT°X°jX 


VK 


VK 


(3.4) 


This means that both thermal boundary layers adjust their thickness and temperature 
drop to yield Ra\ ~ Ra°^ ~ Ra^ — 1000 (e.g., Malkus 1954). The temperature drops 
at both boundaries and the ratio of the thermal boundary layer thicknesses can then be 
derived using Eqs. (3.2-3.3) 


AT* = 


1 


1 -b 773/2 


AT° ~ T^, = 


A? 


1 -b 773/2 yi 


1/4 

Xg 

Tji/a ’ 


where 


9iroy 


(3.5) 


(3.6) 


is the ratio of the gravitational acceleration between the inner and the outer boundaries. 
Figure 4(a) reveals that the marginal stability hypothesis is not fulfilled when different 
radius ratios and gravity prohles are considered. This is particularly obvious for small 
radius ratios where is more than 10 times larger than Ra\. This discrepancy tends 
to vanish when 77 —)■ 1 , when curvature and gravity variations become unimportant. 
As a consequence, there is a significant mismatch between the predicted mean bulk 
temperature from (3.5) and the actual values (figure 4&). Deschamps et al. (2010) also 
reported a similar deviation from (3.5) in their spherical shell models with infinite Prandtl 
number. They suggest that assuming instead Ra‘^/Ra\ ~ 77 ^ might help to improve the 
agreement with the data. This however cannot account for the additional dependence on 
the gravity prohle visible in figure 4. We finally note that Rax < 400 for the database 
of numerical simulations explored here, which suggests that the thermal boundary layers 
are stable in all our simulations. 

Alternatively Wu & Libchaber (1991) and Zhang et al. (1997) proposed that the ther¬ 
mal boundary layers adapt their thicknesses such that the mean hot and cold tempera¬ 
ture fluctuations at mid-depth are equal. Their experiments with Helium indeed revealed 
that the statistical distribution of the temperature at mid-depth was symmetrical. They 
further assumed that the thermal fluctuations in the center can be identified with the 
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Figure 5. (a) Ratio of boundary layer temperature scales (3.7) for various radius ratios and 
gravity profiles. The horizontal dashed line corresponds to the hypothetical identity =6°. 
(b) Temperature drop at the outer boundary layer. The lines correspond to the theoretical 
prediction given in (3.8). 


boundary layer temperature scales 0* = 


and 9° = 


which characterise the 


temperature scale of the thermal boundary layers in a different way than the relative 
temperature drops AT* and AT°. This second hypothesis yields 


6 »* = 9° 


, 3 ’ 


cxpiXip exgoXrp 

and the corresponding temperature drops and boundary layer thicknesses ratio 

1 


AT* = 


1 + rf Xg 


AT° = T^, = 


2 1/3 

r Xg 

1 I 2 

1 + 1? Xg 


\i ' 

/\r^ 


(3.7) 


(3.8) 


Figure 5(a) shows 0°/0* for different radius ratios and gravity profiles, while figure 5(&) 
shows a comparison between the predicted mean bulk temperature and the actual values. 
Besides the cases with g = (ro/r)^ which are in relatively good agreement with the 
predicted scalings, the identity of the boundary layer temperature scales is in general 
not fulfilled for the other gravity profiles. The actual mean bulk temperature is thus 
poorly described by (3.8). We note that previous findings by Ahlers et al. (2006) already 
reported that the theory by Wu & Libchaber’s does also not hold when the transport 
properties depend on temperature (i.e. non-Oberbeck-Boussinesq convection). 


3.3. Conservation of the average plume density in spherical shells 
As demonstrated in the previous section, none of the hypotheses classically employed ac¬ 
curately account for the temperature drops and the boundary layer asymmetry observed 
in spherical shells. We must therefore find a dynamical quantity that could be possibly 
identified between the two boundary layers. 

Figure 6 shows visualisations of the thermal boundary layers for three selected numer¬ 
ical models with different radius ratios and gravity profiles. The isocontours displayed in 
panels (a-c) reveal the intricate plume structure. Long and thin sheet-like structures form 
the main network of plumes. During their migration along the spherical surfaces, these 
sheet-like plumes can collide and convolute with each other to give rise to mushroom- 
type plumes (see Zhou & Xia 20106; Chilla & Schumacher 2012). During this morpho¬ 
logical evolution, mushroom-type plumes acquire a strong radial vorticity component. 
These mushroom-type plumes are particularly visible at the connection points of the 












Figure 6. {a-c) Isosurfaces and equatorial cut of the temperature for three numerical models: 
hot at the inner thermal boundary layer T(r = ri+Ay) in red, cold at the enter thermal boundary 
layer T{r = ro — Ay) in blue, {d-f) Meridional cuts and surfaces of the temperature fluctuations 
T'. The inner (outer) surface corresponds to the location of the inner (outer) thermal boundary 
layers. Color levels range from -0.2 (black) to 0.2 (white). Panels (a) and (d) correspond to 
a model with Ra = 3 x 10®, 77 = 0.3 and g = (ro/r)^- Panels (6) and (e) correspond to a 
model with Ra — 10*, rj = 0.6 and g = (ro/r)^. Panels (c) and (/) correspond to a model with 
Ra = 4 X 10^, g = 0.8 and g = r/ to- 
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T' = a - T' = a/4 . \7„u = 0 


Figure 7. (a) Temperature fluctuation at the inner thermal boundary layer T'{r = ri + Ay) 
displayed in a Hammer projection for a case with Ra = 10®, g = 0.6, g = (ro/r)'^. (b) Cor¬ 
responding binarised extraction of the plumes boundaries using (3.11) and T' < a 12 to define 
the inter-plume area. Zoomed-in contour of the temperature fluctuation T' (c), the horizontal 
divergence V/r • u (d) and the thermal dissipation rate ey (e). The three black contour lines in 
panels (c-e) correspond to several criteria to extract the plume boundaries (3.12). 


sheet plumes network at the inner thermal boundary layer (red isosurface in figure 6a-c). 
Figure 6{d-f) shows the corresponding equatorial and radial cuts of the temperature 
fluctuation T' = T — 'd. These panels further highlight the plume asymmetry between 
the inner and the outer thermal boundary layers. For instance, the case with ry = 0.3 and 
9 — {fo/f)^ (top panels) features an outer boundary layer approximately 4.5 times thicker 
than the inner one. Accordingly, the mushroom-like plumes that depart from the outer 
boundary layer are significantly thicker than the ones emitted from the inner boundary. 
This discrepancy tends to vanish in the thin shell case (ry = 0.8, bottom panels) in which 
curvature and gravity variations play a less significant role (Ay/Ay ~ 1.02 in that case). 

Puthenveettil & Arakeri (2005) and Zhou & Xia (20106) performed statistical analy¬ 
sis of the geometrical properties of thermal plumes in experimental RB convection. By 
tracking a large number of plumes, their analysis revealed that both the plume separation 
and the width of the sheet-like plumes follow a log-normal probability density function 
(PDF). 

To further assess how the average plume properties of the inner and outer thermal 
boundary layers compare with each other in spherical geometry, we adopt a simpler 
strategy by only focussing on the statistics of the plume density. The plume density per 
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surface unit at a given radius r is expressed by 

N 


(3.9) 


where N is the number of plumes, approximated here by the ratio of the spherical surface 
area to the mean inter-plume area S: 


iV~ 




(3.10) 


This inter-plume area S can be further related to the average plume separation i via 

5- (7r/4)P. 

An accurate evaluation of the inter-plume area for each thermal boundary layer how¬ 
ever requires to separate the plumes from the background fluid. Most of the criteria 
employed to determine the location of the plume boundaries are based on thresholds of 
certain physical quantities (see Shishkina & Wagner 2008, for a review of the different 
plume extraction techniques). This encompasses threshold values on the temperature 
fluctuations T' (Zhou & Xia 2002), on the vertical velocity Ur (Ching et al. 2004) or on 
the thermal dissipation rate ex (Shishkina & Wagner 2005). The choice of the threshold 
value however remains an open problem. Alternatively, Vipin & Puthenveettil (2013) 
show that the sign of the horizontal divergence of the velocity V/f • it might provide a 
simple and threshold-free criterion to separate the plumes from the background fluid 


Vh ■ u = 


1 


d 


r sin 6 d6 


(sin Oug) + 


duA 


r sin 6 d(j) 


J.2 


Ur 


Fluid regions with V// • it < 0 indeed correspond to local regions of positive vertical accel¬ 
eration, expected inside the plumes, while the fluid regions with Vjj • m > 0 characterise 
the inter-plume area. 

To analyse the statistics of S, we thus consider here several criteria based either on a 
threshold value of the temperature fluctuations or on the sign of the horizontal divergence. 
This means that a given inter-plume area at the inner (outer) thermal boundary layer is 
either defined as an enclosed region surrounded by hot (cold) sheet-like plumes carrying 
a temperature perturbation \T'\ > t; or by an enclosed region with Vjj - it > 0. To further 
estimate the possible impact of the chosen threshold value on S, we vary t between a/A 
and cr. This yields 

S{r)=r^ <j) sin0ddd()), (3-11) 

where the physical criterion 7) (To) to extract the plume boundaries at the inner (outer) 
boundary layer is given by 


T = 

To = 


T'(r\, te [a(rl),a(r\/2),a(r\/A)], 

• M > 0, 

T'(rl, e,(j))>t, t€ [(7(r\), a(r\/2), a(r\/A)], 

V_f/ ■ M > 0, 


(3.12) 


where r\ = ri + (r^ = Tq — A^) for the inner (outer) thermal boundary layer. 

Figure 7 shows an example of such a characterisation procedure for the inner thermal 
boundary layer of a numerical model with Ra = 10®, rj = 0.6, g = (rojr)'^. Panel (b) 
illustrates a plume extraction process when using \T\ > a 12 to determine the location 
of the plumes: the black area correspond to the inter-plume spacing while the white area 
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Figure 8. Probability density functions (PDFs) of the dimensionless inter-plume area 5 at the 
inner (a) and at the outer (6) thermal boundary layers using different criteria to extract the 
plumes (3.12) for a model with Ra = 4 X 10^ rj = 0.8 and g = r/r^. 


correspond to the complementary plume network location. The fainter emerging sheet¬ 
like plumes are filtered out and only the remaining “skeleton” of the plume network 
is selected by this extraction process. The choice of tT/2 is however arbitrary and can 
influence the evaluation of the number of plumes. The insets displayed in panels (c-e) 
illustrate the sensitivity of the plume extraction process on the criterion employed to 
detect the plumes. In particular, using the threshold based on the largest temperature 
fluctuations \T'\ > cr can lead to the fragmentation of the detected plume lanes into 
several isolated smaller regions. As a consequence, several neighbouring inter-plume areas 
can possibly be artificially connected when using this criterion. In contrast, using the sign 
of the horizontal divergence to estimate the plumes location yields much broader sheet¬ 
like plumes. As visible on panel (e), the plume boundaries frequently correspond to local 
maxima of the thermal dissipation rate ct (Shishkina & Wagner 2008). 

For each criterion given in (3.12), we then calculate the area of each bounded black 
surface visible in figure 7(6) to construct the statistical distribution of the inter-plume 
area for both thermal boundary layers. Figure 8 compares the resulting PDFs obtained 
by combining several snapshots for a numerical model with Ra = 4 x 10^, rj = 0.8 and 
g = rlro. Besides the criterion \T'\ > a which yields PDFs that are slightly shifted 
towards smaller inter-plume spacing areas, the statistical distributions are found to be 
relatively insensitive to the detection criterion (3.12). We therefore restrict the following 
comparison to the criterion \T'\ > a 12 only. 

Figure 9 shows the PDFs for the three numerical models of figure 6. For the two cases 
with r\ = 0.6 and r\ = 0.8 (panels 6-c), the statistical distributions for both thermal 
boundary layers nearly overlap. This means that the inter-plume area is similar at both 
spherical shell surfaces. In contrast, for the case with r] = 0.3 (panel a), the two PDFs 
are offset relative to each other. However, the peaks of the distributions remain relatively 
close, meaning that once again the inner and the outer thermal boundary layers share a 
similar average inter-plume area. Puthenveettil & Arakeri (2005) and Zhou & Xia (20106) 
demonstrated that the thermal plume statistics in turbulent RB convection follow a 
log-normal distribution (see also Shishkina & Wagner 2008; Puthenveettil et al. 2011). 
The large number of plumes in the cases with rj = 0.6 and rj = 0.8 (hgure 66-c) would 
allow a characterisation of the nature of the statistical distributions. However, this would 
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Figure 9. PDF of the dimensionless inter-plume area S at the outer (orange upward triangles) 
and at the inner (blue downward triangles) thermal boundary layers using \T'\ > a 12 to extract 
plumes. Panel (a) corresponds to a numerical model with Ra = 3 x 10®, tj = 0.3 and g = (ro/r)®. 
Panel (b) corresponds to a model with Ra = 10®, g — 0.6 and g = (ro/r)^. Panel (c) corresponds 
to a model with Ra = 4 x 10^, g = 0.8 and g = r/ro- The two vertical lines correspond to the 
predicted values of the mean inter-plume area 5 derived from (3.15) for both thermal boundary 
layers. 



X 



Figure 10. Schematic showing two adjacent plumes separated by a distance The thick black 
arrows indicate the direction of merging of the two plumes. 




be much more difficult in the g = 0.3 case (figure 6a) in which the plume density is 
significantly weaker. As a consequence, no further attempt has been made to characterise 
the exact nature of the PDFs visible in figure 9, although the universality of the log¬ 
normal statistics reported by Puthenveettil & Arakeri (2005) and Zhou & Xia (20106) 
likely indicates that the same statistical distribution should hold here too. 

The inter-plume area statistics therefore reveals that the inner and the outer thermal 
boundary layers exhibit a similar average plume density, independently of the spherical 
shell geometry and the gravity profile. Assuming ~ would allow us to close the 
system of equations (3.2-3.3) and thus finally estimate AT*, AT° and A^/A^. This how¬ 
ever requires us to determine an analytical expression of the average inter-plume area 
S or equivalently of the mean plume separation t that depends on the boundary layer 
thickness and the temperature drop. 
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Using the boundary layer equations for natural convection (Rotem & Claassen 1969), 
Puthenveettil et al. (2011) demonstrated that the thermal boundary layer thickness fol¬ 
lows 


\ 2,0 
y\rp 


{x) 


(Rar)i/5’ 


(3.13) 


where x is the distance along the horizontal direction and Ratf = agAT'^'°x^ /vk is a 
Rayleigh number based on the lengthscale x and on the boundary layer temperature 
jumps As shown on figure 10, using x = (Puthenveettil & Arakeri 2005; 

Puthenveettil et al. 2011) then allows to establish the following relation for the average 
plume spacing 


which yields 


Xt 1 
^ RaY^' 




agiAT^Xlp 


'agoAT°X^^ 


(3.14) 


(3.15) 


for both thermal boundary layers. We note that an equivalent expression for the average 
plume spacing can be derived from a simple mechanical description of the equilibrium 
between production and coalescence of plumes in each boundary layer (see Parmentier 
& Sotin 2000; King et al. 2013). 

Equation (3.15) is however expected to be only valid at the scaling level. The vertical 
lines in figure 9 therefore correspond to the estimated average inter-plume area for both 
thermal boundary layers using (3.15) and Si^o — The predicted average inter¬ 

plume area is in good agreement with the peaks of the statistical distributions for the 
three cases discussed here. The expression (3.15) therefore provides a reasonable estimate 
of the average plume separation (Puthenveettil & Arakeri 2005; Puthenveettil et al. 2011; 
Gunasegarane & Puthenveettil 2014). The comparable observed plume density at both 
thermal boundary layers thus yields 



ag.AT^Xij.^ _ agoAT°X^^^ 

un UK 


(3.16) 


Using Eqs. (3.2-3.3) then allows us to finally estimate the temperature jumps and the 
ratio of the thermal boundary layer thicknesses in our dimensionless units: 


AT* = 


1 

1 -h r?5/3 ’ 


AT° = Trr 


„5/3..^,1/6 xo 1/6 

V Xg ^ Xg 

1 _(_ jy5/3 7^1/3 


(3.17) 


Figure 11 shows the ratios X^/X\, and the temperature jumps AT* and AT°. 

In contrast to the previous criteria, either coming from the marginal stability of the 
boundary layer (3.5, figure 4) or from the identity of the temperature fluctuations at 
mid-shell (3.17, figure 5), the ratio of the average plume separation £o/G now falls much 
closer to the unity line. Some deviations are nevertheless still visible for spherical shells 
with T] < 0.4 and g = r/vo (orange circles). The comparable average plume density 
between both boundary layers allows us to accurately predict the asymmetry of the 
thermal boundary layers A^/A^ and the corresponding temperature drops for the vast 
majority of the numerical cases explored here (solid lines in panels b-d). 

As we consider a fluid with Pr = 1, the viscous boundary layers should show a com- 
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Figure 11. (a) Ratio of the thermal plume separation estimated by (3.15) for various radius 
ratios and gravity profiles. The horizontal dashed line corresponds to the identity of the aver¬ 
age plume separation between both thermal boundary layers, i.e. t = . (b) Ratio of thermal 

boundary layer thicknesses, (c) Temperature drop at the inner boundary layer, (d) Tempera¬ 
ture drop at the outer boundary layer. The lines in panels (b-d) correspond to the theoretical 
prediction given in (3.17). 



Figure 12. Ratio of the viscous boundary layer thicknesses for various aspect ratios and 
gravity proHles. The lines correspond to the theoretical prediction given in (3.18). 


parable degree of asymmetry as the thermal boundary layers. (3.17) thus implies 


^ _ Xg^ 
xIj " “ 77 '/^ ■ 


(3.18) 


Figure 12 shows the ratio of the viscous boundary layer thicknesses for the different 
setups explored in this study. The observed asymmetry between the two spherical shell 
surfaces is in a good agreement with (3.18) (solid black lines). 
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Figure 13. (a) Thermal boundary layer thicknesses at the outer boundary (Ay) and at the 
inner boundary (A^) as a function of the Nusselt number for the cases of Table 1. The two lines 
correspond to the theoretical predictions from (3.19). (6) Normalised boundary layer thicknesses 
as a function of the Nusselt number for different radius ratios and gravity profiles. For the sake 
of clarity, the outer boundary layer is only displayed for the cases with g = (rojr'f'. The solid 
line corresponds to At = 0.5 Nu~^ (3.20). 


3.4. Thermal boundary layer scalings 

Using (3.17) and the definition of the Nusselt number (2.6), we can derive the following 
scaling relations for the thermal boundary layer thicknesses: 


A^ = 


1 + 


1 ^ 

1 + rf!'^ xj^ 


(3.19) 


Figure 13(a) demonstrates that the boundary layer thicknesses for the numerical sim¬ 
ulations of Table 1 {g = (ro/r)^ and rj = 0.6) are indeed in close agreement with the 
theoretical predictions. To further check this scaling for other spherical shell configura¬ 
tions, we introduce the following normalisation of the thermal boundary layer thicknesses 


^ 2 r, 2 ^ 2/3 


This allows us to derive a unified scaling that does not depend on the choice of the gravity 
profile or on the spherical shell geometry 


\T = Xi=X° = 


2Nu 


(3.20) 


Figure 13(&) shows this normalised boundary layer thickness for the different spherical 
shell configurations explored here. Despite the variety of the physical setups, the nor¬ 
malised boundary layer thicknesses are in good agreement with the predicted behaviour. 
A closer inspection of figure 11(&) and Table 1 nevertheless reveals a remaining weak 
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Figure 14. Ratio of the thermal boundary layer thicknesses A^/Ay as a function of the Rayleigh 
number for the numerical models of Table 1 with rj = 0.6 and g = (ro/r)^. The horizontal dashed 
line corresponds to the predicted ratio A^/A^ given in (3.17). 

dependence of the ratio A^/A^ on the Rayleigh number. This is illustrated in figure 14 
which shows A^/A^ as a function of Ra for the ij = 0.6, g = {vo/r)^ cases of Table 1. 
Although some complex variations are visible, the first-order trend is a very slow increase 
of X^/X^rp with Ra (10% increase over five decades of Ra). No evidence of saturation is 
however visible and further deviations from the predicted ratio (horizontal dashed line) 
might therefore be expected at larger Ra. This variation might cast some doubts on the 
validity of the previous derivation for higher Rayleigh numbers. This may imply that 
either the plume separation is not conserved at higher Ra; or that the estimate of the 
average plume spacing is too simplistic to capture the detailed plume physics in turbulent 
convection. 


4. Boundary layer analysis 

The Grossmann Lohse (GL) theory relies on the assumption that the viscous and the 
thermal boundary layers are not yet turbulent. This is motivated by the observation of 
small boundary layer Reynolds numbers Rcs = Re X/d < 200 in experimental convection 
up to Ra ~ 10^^, which remain well below the expected transition to fully turbulent 
boundary layers (expected at Reg ~ 420, see Ahlers et al. 2009). The boundary layer 
flow is therefore likely laminar and follows the Prandtl-Blasius (PB) laminar boundary 
layer theory (Prandtl 1905; Blasius 1908; Schlichting & Gersten 2000). The PB theory 
assumes a balance between the viscous forces, important in the boundary layers, and 
inertia which dominates in the bulk of the fluid. For the numerical models with unity 
Prandtl number, this directly implies that the boundary layer thicknesses are inversely 
proportional to the square-root of Re 

Xjj ~ Xt ^ Re . (4-1) 

Figure 15(a-6) shows a test of this theoretical scaling for the two different definitions 
of the viscous boundary layer introduced in § 3.1. Confirming previous findings by Breuer 
et al. (2004), the commonly-employed definition based on the location of the horizontal 
velocity maxima yields values that significantly differ from the theoretical prediction 
(4.1). The least-square fit to the data for the cases with Re > 250 indeed gives values 
relatively close to A( 7 ,m ~ an exponent already reported in the experiments by 

Lam et al. (2002) and in the numerical models in cartesian geometry by Breuer et al. 
(2004) and King et al. (2013). In addition, A^.m is always significantly larger than Xt, 
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Figure 15. (a) Viscous boundary layer thicknesses at the outer boundary {\\j^rn) at the 
inner boundary {Xu,m) as a function of the Reynolds number for the cases of Table 1 (17 = 0.6, 
g = (vo/r)'^). (b) Same for the other definition of the viscous boundary layer, i.e. AJj and A^. 
(c) Thermal boundary layer thicknesses at the outer boundary (Ay) and at the inner boundary 
(Ay) as a function of Re. The black lines in panels (a, b and c) correspond to the least-square 
fit to the data for the numerical models with Ra > 10® (i.e. Re > 250). 


which is at odds with the expectation At — Ac/ when Pr = 1 (see Table 1 for detailed 
values). 

Adopting the intersection of the two tangents to define the viscous boundary layers 
(figure 156) leads to exponents much closer to the predicted value of 1/2 in the high-i?e 
regime. The viscous boundary layer thicknesses obtained with this second definition are 
in addition found to be relatively close to the thermal boundary layer thicknesses in the 
high Reynolds regime, i.e. Xu ~ At- Therefore, both the expected scaling of Xu with 
Re and the similarities between thermal and viscous boundary layer thicknesses strongly 
suggest that the tangent-intersection method is a more appropriate way to estimate the 
actual viscous boundary layer. We therefore focus on this definition in the following. 
For low Reynolds numbers {Re < 200), however, the viscous boundary layer thicknesses 
deviate from the PB theory and follow a shallower slope around This deviation 

implies a possible inaccurate description of the low Ra cases by the GL theory (see 
below). 

Figure 15(c) shows that the corresponding scaling of the thermal boundary layer with 
Re follows a similar trend as the viscous boundary layers. The best fit to the data for the 
cases with Re > 250 yields exponents moderately larger than the theoretical prediction 
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Figure 16. {a-h) Radial profiles of the time and horizontally-averaged temperature il and 

horizontal velocity Rch- (c) 0 as a function of ^t- (d) Uh as a function of The solid lines 
in panels (c-d) corresponds to the Prandtl-Blasius solution. The inset in panel (c) shows 0 in 
double-logarithmic scale. For the sake of clarity, the outer boundary layer is only displayed for 
one single case in the panels (c-d) {Ra — 10®, 77 = 0.6, g = (ro/r)®). 


(4.1), while the low—i?e cases follow a shallower exponent. At this point, we can simply 
speculate that this difference might possibly arise because of the inherent dynamical 
nature of the thermal boundary layers. 

For a meaningful comparison with the boundary layer theory, we define new scaling 
variables for the distance to the spherical shell boundaries. These variables are introduced 
to compensate for the changes in the boundary layer thicknesses that arise when i?a, rj 
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or g are modified. This then allows us to accurately characterise the shape of both the 
temperature and the flow profiles in the boundary layers and to compare them with the 
PB boundary layer profiles. To do so, we introduce the self-similarity variables and 
for both the inner and the outer spherical shell boundaries: 


— 


r — ri 

Xl 


Ct — 


To - r 

X°T ' 


Cu = 


r - n 

Xlr 


Cu = 


To-r 
\° 


'T ''T ''U 

We accordingly define the following rescaled temperatures for both boundaries 


f\r,e,cl),t) = 


T{r,9,(j),t) - -dirm) 
Tbot - 


f°{r,e,cj>,t) = 


'd(r^) - T{r,9,4i,t) 

- Ttop 


(4.2) 


(4.3) 


where r™ = {ri ro)l2 is the mid-shell radius. The rescaled horizontal velocity is simply 
obtained by normalisation with its local maximum for each boundary layer: 


u'‘h{r,9,4),t) 


Uh{r, 9 ,(j),t) 

iaaxi{Reh) ’ 


Uh{r,9,(l),t) 


Uh{r,9,(j),t) 
maxo(i?e;i) ‘ 


(4.4) 


To check the similarity of the profiles, we consider five numerical models with different Ra, 
rj and g. Figure 16(a-6) show the typical mean horizontal velocity and temperature for 
these cases, while figure 16(c-d) show the corresponding time and horizontally-averaged 
normalised quantities: 

0 = , Ub = (ur)^- (4-5) 

As already shown in the previous section, the bulk temperature as well as the bound¬ 
ary layer asymmetry strongly depend on the gravity profile and the radius ratio of the 
spherical shell. Increasing Ra leads to a steepening of the temperature profiles near the 
boundaries accompanied by a shift of the horizontal velocity maxima towards the walls. 
Although d and Ret drastically differ in the five cases considered here, introducing the 
normalised variables 0 and lAh allows to merge all the different configurations into one 
single radial profile. The upward and downward pointing triangles further indicate that 
those profiles are also independent of the choice of the boundary layer (at the inner or 
at the outer boundary). Finally, the solutions remain similar to each other when Ra is 
varied, at least in the interval considered here (i.e. 10® < Ra < 10®). These results are 
in good agreement with the profiles obtained in the numerical simulations by Shishkina 
& Thess (2009) that cover a similar range of Rayleigh numbers in cylindrical cells with 
F = 1. 

Figure 16(c-d) also compares the numerical profiles with those derived from the PB 
boundary layer theory. The time-averaged normalised temperature and velocity profiles 
slightly deviate from the PB profiles, confirming previous 2-D and 3-D numerical models 
by Zhou & Xia (2010a), Shishkina & Thess (2009) and Stevens et al. (2010) This deviation 
can be attributed to the intermittent nature of plumes that permanently detach from the 
boundary layers (Sun et al. 2008; Zhou & Xia 2010a; du Puits et al. 2013). When the 
boundary layer profiles are obtained from a time-averaging procedure at a fixed height 
with respect to the container frame (i.e. and ^jj are time-independent), they can 
actually sample both the bulk and the boundary layer dynamics as the measurement 
position can be at time either inside or outside the boundary layer (for instance during 
a plume emission). To better isolate the boundary layer dynamics, Zhou & Xia (2010a) 
therefore suggested to study the physical properties in a time-dependent frame that 
accounts for the instantaneous boundary layer fluctuations (see also Zhou et al. 2010; 
Stevens et al. 2012; Shishkina et al. 2015). 

We apply this dynamical rescaling method to our numerical models by defining local 
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Figure 17. (a) Time and horizontally-averaged normalised temperature profile in the fixed 

reference frame (dotted line, 4.5) and in the dynamical frame (dashed line, 4.6) for a case with 
Ra = 10®, 77 = 0.6 and g = (uo/r)^. (6) Corresponding horizontal velocity profile in the fixed 
reference frame (dotted line, 4.5) and in the dynamical frame (dashed line, 4.7). The solid lines 
in both panels correspond to the Prandtl-Blasius solution. 


and instantaneous boundary layer thicknesses 




' O ' 

X^(0, (j), t) 






As the inner and the outer boundary layers exhibit the same behaviour (figure 16), we 
restrict the following discussion to the outer boundary layer. Following Zhou et al. (2010) 
and Shi et al. (2012), the horizontal velocity and temperature profiles are given by 


i^u) = {^h (r, = r, 




S ’ 


(4.6) 


0*(^J) = (r,0,<^,t|r = ro-^*A§.(d,((.,t)))^. (4.7) 

Practically, this dynamical rescaling strategy has been achieved by measuring the local 
and instantaneous boundary layer for each grid coordinates for several snapshots. 

Following Zhou & Xia (2010a) and Stevens et al. (2012), the temperature profiles have 
been further normalised to some position outside the thermal boundary layer (here = 5 
or = 5) to ease the comparison with the classical definition of the boundary layer in 
the fixed reference frame (4.5). Figure 17 shows an example of this dynamical rescaling 
method applied to a case with Ra = 10®, t] = 0.6 and g = (ro/r)^. The temperature and 
horizontal velocity profiles in the spatially and temporally varying local frame are now 
in much closer agreement with the PB laminar profiles than those obtained in the fixed 
reference frame. 

Because of the numerical cost of the whole procedure, the dynamical rescaling has only 
been tested on a limited number of cases. Applying the same method to the numerical 
model with Ra = 10® {rj = 0.6, g = (ro/r)^) yields nearly indistinguishable profiles from 
those displayed in figure 17. This further indicates that boundary layers in spherical shells 
are laminar in the Ra range explored here and can be well described by the PB theory, 
provided boundary layers are analysed in a time-dependent frame, which fluctuates with 
the local and instantaneous boundary layer thicknesses. 
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Figure 18. Measured contributions of the boundary layer (open symbols) and the fluid bulk 
(filled symbols) to the total viscous dissipation rate £(/ (orange squares) and the total thermal 
dissipation rate et (blue circles). 

5. Dissipation analysis 

5.1. Bulk and boundary layer contributions to viscous and thermal dissipation rates 
The prerequisite of a laminar boundary layer seems fulfilled in our numerical models and 
we can thus try to apply the GL formalism to our dataset. The idea of the GL theory is 
to separate the viscous and thermal dissipation rates into two contributions, one coming 
from the fluid bulk (indicated by the superscript bu in the following) and one coming 
from the boundary layers (bl), such that 

£7' = £7^+6^, £j/ = £^+£y, ( 5 . 1 ) 

where the contributions from the bulk are defined by 



and the boundary layer contributions are given by 



The RB flows are then classified in the Ra — Pr parameter space according to the 
dominant contributions to the viscous and thermal dissipation rates. This defines four 
regimes depending on Ra and Pr: regime 1 when eu ~ and ct — regime II when 
£[/ ~ £^ and £t — £^1 regime III when eu — e^ and ct — £^; and finally regime IV 
when eu ~ £^ and ct e^. 

For a unity Prandtl number, the GL theory predicts that the flows should be dominated 
by dissipations in the boundary layer regions at low Ra (regime I). A transition to another 
regime where dissipations in the fluid bulk dominate (regime IV) is expected to happen 
roughly around Ra ~ 10® — 10^° (Grossmann & Lohse 2000; Ahlers et al. 2009; Stevens 
et al. 2013). 

Figure 18 shows the relative contributions of the bulk and boundary layers to the 
viscous and thermal dissipation rates. The viscous dissipation eu is always dominated 
by the bulk contribution: starting from roughly 60% at Ra = 10"^, it nearly reaches 
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90% at Ra = 10®. In contrast, the thermal dissipation rate is always dominated by the 
boundary layer regions, such that slowly increases from 10% at Ra = 10^ to 30% at 
Ra = 10®. According to the GL classification, all our numerical simulations thus belong 
to the regime II of the Ra — Pr parameter space, in which and > e^. This 

seems at odds with the GL theory, which predicts that our DNS should be located either 
in the regime I or in the regime IV of the parameter space for the range of Ra explored 
here (10^ < Ra < 10®). 

The dominance of the boundary layer contribution in the thermal dissipation rate was 
already reported by Verzicco (2003) for the same range of Ra values. This phenomenon 
may be attributed to the dynamical nature of the plumes which permanently detach 
from the boundary layers and penetrate in the bulk of the fluid. These thermal plumes 
have the same typical size as the boundary layer thickness and can thus be thought 
as “detached boundary layers”. Grossmann & Lohse (2004) have therefore suggested to 
modify their scaling theory to incorporate these detached boundary layers in the thermal 
dissipation rate. They propose to decompose ex into one contribution coming from the 
plumes (cy) and one coming from the turbulent background (e^) 

eT = er+ef', (5.2) 

Such a decomposition is however extremely difficult to conduct in spherical shells in which 
the very large aspect ratio of the convective layer yields a complex and time-dependent 
multi-cellular large scale circulation (LSG) pattern (see for instance Bailon-Guba et al. 
2010, for the influence of large T on the LSG). In the following, we therefore first keep 
the initial decomposition (5.1) before coming back to the inherent problem of accurately 
separating the different contributions to the dissipation rate in § 5.3. 


5.2. Individual scaling laws for the dissipation rates 

Based on the hypothesis of homogeneous and isotropic turbulence, the GL theory assumes 
that the thermal and viscous dissipation rates in the bulk of the fluid scale like 

e^“ ~ i?e^ ~ Re, (5.3) 


in our dimensionless units. Figure 19 shows the bulk dissipation rates as a function of 
Re for the numerical models of Table 1. The best fit to the data (solid lines) for the 
cases with Ra > 10^ yields e’ff ^ i?e®'^® and ^ only roughly similar to the 

prediction (5.3). These deviations from the theoretical exponents are further confirmed 
by the compensated scalings i?e“® and Re~^ shown in panels (b) and (d), which 
show a coherent remaining dependence on Re. Even at high Reynolds numbers, there is 
no evidence of convergence towards the exact expected scalings from the GL theory. This 
is particularly obvious for which remains in close agreement with ^ for the 

whole range of Re values explored here (solid line in figure 19d). The dependence of 
on Re shows a gradual steepening of the slope when Re increases, which implies that 
e\j‘{Re) cannot be accurately represented by a simple power law. Similar deviations from 
(5.3) have already been reported in the Taylor-Gouette flow experiments by Lathrop et al. 
(1992) and in the numerical simulations of RB homogeneous turbulence by Galzavarini 
et al. (2005). 

A similar procedure can be applied to the dissipation rates in the boundary layers. 
In spherical shells, the volume fraction occupied by the boundary layers can be approxi¬ 
mated by the following Taylor expansion to the first order 
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Figure 19. (a) Viscous dissipation in the bulk of the fluid as a function of Re. {b) Corresponding 
compensated scaling, (c) Thermal dissipation in the bulk of the fluid as a function of Re. 
(d) Corresponding compensated el^ scaling. The solid black lines in the four panels correspond 
to the least-square fit to the data for the numerical models with Ra > 10®. 


in the limit of thin boundary layers A <C ro — Vi. The viscous dissipation rate in the 
boundary layers can then be estimated by 

Tj2 

bl rms £ f\ \ 

~ -.2 ~ “T , 

A (7 -\u 

in our dimensionless units. As demonstrated in § 4, the boundary layers are laminar 
and are in reasonable agreement with the PB boundary layer theory. This implies that 
\u ^ and thus yields 

(5.4) 

Figure 20 shows the viscous dissipation rate in the boundary layers as a function of Re 
for the numerical models of Table 1. The least-square fit to the data (solid lines) yields 
Cy ^ in close agreement with the expected theoretical exponent. The compensated 

scaling displayed in panel (6) reveals a remaining weak secondary dependence of on 
Re, which is not accurately captured by the power law. The local slope of e’lj{Re) initially 
decreases with Re (when Re < 200) before slowly increasing at higher Reynolds numbers 
{Re > 200). This behaviour suggests that, although simple power laws are at first glance 
in very good agreement with the GL theory, they may not account for the detailed 
variations of e^{Re). 

The boundary layer contribution to the thermal dissipation rate is estimated in a 
similar way: 

M ^ AT2 1 

^ ^ ^ Xt' 
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Figure 20. (a) Viscous dissipation in the boundary layers as a function of Re. {b) Corresponding 
compensated scaling. The solid black lines in the four panels correspond to the least-square 
fit to the data for the numerical models with Ra > 10®. 



Figure 21. (a) Thermal dissipation in the boundary layers as a function of Nu. (b) Correspond¬ 
ing compensated scaling, (c) Thermal dissipation in the boundary layers as a function of Re. 
(d) Corresponding compensated scaling. The solid black lines in the four panels correspond 
to the least-square fit to the data for the numerical models with Ra > 10®. 

which yields 

~ Nu. (5.5) 

The laminar nature of the boundary layers also implies Ay ^ and thus 

Ct ~ (5.6) 

Figure 21 shows as a function of Nu and Re for the cases of Table 1. The least-square 
fits yield ^ and ^ (solid lines), close to the expected exponents. 

However, the compensated scalings displayed in panels (6) and (d) reveal that the linear 
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fit to e^{Nu) remains in good agreement with the data, while e^{Re) is not accurately 
described by such a simple fit. The solutions increasingly deviate from the power law at 
high Reynolds numbers with the local slope of e^{Re) that gradually steepens with Re. 


5.3. Individual versus global scalings 

Despite the overall fair agreement with the GL predictions, a close inspection of the de¬ 
pendence of the four dissipation rates on the Reynolds number reveals some remaining 
dependence on Re, which cannot be perfectly described by simple power laws. This is 
particularly obvious in the boundary layer contributions e’^(Re) and e^{Re). In addition, 
both thermal dissipation rates deviate stronger from the theoretical exponents than their 
viscous counterparts (see also Verzicco 2003). One obvious problem is the inherent dif¬ 
ficult separation of bulk and boundary layer contributions already discussed above. The 
dynamical plumes constantly departing from the boundary layers obviously complicate 
matters. To check whether the general idea of a boundary layer and a bulk contribu¬ 
tion that both scale with the predicted exponents is at least compatible with the total 
dissipation rates, we directly fit 

?)) = -f = a Re^ -b b , 

ef = = cRe-\- dRe^^^ . 

This leaves only the four prefactors (a, 6, c, d) as free fitting parameters and yields 

^ = 0.248Re3-b 7.084 

Z. Z. (5.7) 

?^ = e^“-be^= 0.004 Re-b 0.453 Re^/^ 


We then compare this direct least-square fit of the total dissipation rates to the sum of 
the individual scalings obtained in the previous section 


= e^“ -b = 1.756 Re^-'^® -b 2.197 Re^ '^^ , 


_ ^bu 

CT — Crp 


= 0.038 Re°-^ -b 0.268 Re° ®^ , 


(5.8) 


The accuracy of the two scalings (5.7) and (5.8) are compared in figure 22, which shows 
the total viscous and thermal dissipation rates as a function of Re for the numerical 
models of Table 1. While the two scalings are nearly indistinguishable on the left panels 
(a) and (c), the corresponding normalised scalings shown in panels (6) and (d) reveal 
some important differences. The scalings based on the sum of the power laws derived in 
the last section (5.8) are in relatively poor agreement with the data (5 — 10% error for 
Re > 10^) with no obvious asymptotic behaviour. On the other hand, the global scalings 
(5.7) fall much closer to the actual values for the range 10^ < Re < 10"*^ and approach an 
asymptote for Re > 10^. The deviations observed for the highest Re cases have probably 
a numerical origin; the averaging timespan used to estimate the dissipation rates are 
likely too short in the most demanding cases to perfectly average out all the fluctuations. 

The total thermal and viscous dissipation rates in our spherical shell simulations are 
thus better described by the sum of two power laws that follow the GL theory than by 
the sum of the asymptotic laws derived from the individual contributions, which suffers 
from an unclear separation of the boundary layer and bulk dynamics. 

This result also sheds a new light on the placement of our numerical simulations in the 
GL regime diagram (figure 18). Equation (5.7) directly provides the estimated relative 
contributions of the bulk and boundary layers to the viscous and thermal dissipation 
rates. Figure 23 shows these different contributions for the numerical models of Table 1 
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Figure 22. (a) tu as a function of Re. (b) eu normalised by the predictions coming from (5.7) 
(orange triangles) and (5.8) (blue circles) as a function of Re. (c) ex as a function of Re. (d) tr 
normalised by the predictions coming from (5.7) (orange triangles) and (5.8) (blue circles) as a 
function of Re. The dashed black lines in panels (c-d) correspond to the equality between the 
asymptotic scalings and the data. 



Figure 23. Estimated relative contributions of the boundary layer (open symbols) and the 
fluid bulk (filled symbols) to eu (orange squares) and er (blue circles) using (5.7). 


and reveals a completely different balance than in figure 18. At low Rayleigh numbers, 
the estimated boundary layer contributions now dominate both the viscous and thermal 
dissipation rates. The viscous dissipation rate in the fluid bulk gradually increases with 
Ra and dominates beyond Ra > 10^. The bulk contribution to the thermal dissipation 
rate exhibits a similar trend, gradually increasing from roughly 5% at Ra = 10^ to 
more than 40% at Ra = 10®. While the thermal dissipation rate in the fluid bulk never 
dominates in the regime explored here, our scaling predict that it will do so for Ra > 
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3 X 10®. These values would then locate the numerical models with Ra < 10^ in the GL 
regime I of the Ra — Pr parameter space. The cases with 10^ < Ra < 3 x 10® would then 
belong to regime II. The transition to regime IV where the bulk contributions dominate 
the dissipation rates would then happen around Ra ~ 3 x 10®. 

While it is not clear that the contributions inferred from the total dissipation rates 
actually reflect the exact bulk and boundary layer contributions to ejj and ct, this separa¬ 
tion nevertheless allows us to reconcile the classihcation of our spherical shell convection 
models in the Ra — Pr parameter space with the prediction of the GL theory for a fluid 
with Pr = 1. Future work that will better characterise and separate the bulk and bound¬ 
ary layer dynamics might help to reconcile the individual scalings (figure 18) with the 
global ones (hgure 23), especially at low Ra. In particular, considering dissipation layers 
as defined by Petschel et al. (2013) instead of classical boundary layers might possibly 
help to better separate the bulk and boundary layer contributions. 


6. Nusselt and Reynolds numbers scalings 


Figure 24 shows Nu and Re as a function of Ra for the cases of Table 1. A simple best fit 
to the data for the cases with Ra > 10® yields Nu ~ and Re ^ iia® "^’^®, relatively 

close to Nu ^ Ra?/"^ and Re ^ Ra^^^. While reducing the scaling behaviours of Nu and 
Re to such simple power laws is a common practice in studies of convection in spherical 
shells with infinite Prandtl number (e.g. Wolstencroft et al. 2009; Deschamps et al. 2010), 
this description might be too simplistic to account for the complex dependence of Vu 
and Re upon Ra. To illustrate this issue, the panels (c) and (d) of figure 24 show 
the compensated scalings of Nu and Re. The power laws fail to capture the complex 
behaviour of Nu{Ra) and Re{Ra) and show an increasing deviation from the data at 
high Ra. For instance, the 0.289 scaling exponent obtained for the Nusselt number is too 
steep for Ra < 10^ and too shallow for higher Rayleigh numbers. 

The GL theory predicts a gradual change of the slopes of Nu{Ra) and Re{Ra) since 
the flows cross different dynamical regimes when Ra increases. Using the asymptotic laws 
obtained for the different contributions to the dissipation rates (5.8) and the dissipation 
relations (2.9) and (2.10), we can derive the following equations that relate Nu and Re 
to Ra-. 


££/ = 


Ct = 


1 -\- rj rj^ Pr^ 


^{Nu- 1) =1.756i?e®-^®-b2.197i?e 


„2.52 


3r] 


■ Nu 


=0.038 -t 0.268 


( 6 . 1 ) 


1 -\-r] -\-r]^ 

This system of equations can be numerically integrated to derive the scaling laws for Nu 
and Re. Figure 24(c-d) illustrates the comparison between these integrated values and 
the actual data, while hgure 25 shows the local effective exponents aeff and (3eg of the 
Nu{Ra) and Re{Ra) laws as a function of Ra: 


d In Nu cl In Re 

91ni?a’ d\nRa ' 


( 6 . 2 ) 


While Re{Ra) is nicely described by the solution of (6.1), some persistent deviations in 
the Nu{Ra) scaling are noticeable. In particular, aeQ{Ra) increases much faster than 
expected from the scaling law (dashed lines): aeff(10®) ~ 0.32 while the predicted slope 
remains close to 0.285. The difficulties to accurately separate the bulk and boundary layer 
dynamics when deriving the scaling laws for the different contributions to the dissipation 
rates are once again likely responsible of this misht. 

As demonstrated in the previous section, replacing the sum of the individual dissipation 
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Figure 24. (a) Nusselt number versus Rayleigh number, (b) Reynolds number versus Rayleigh 
number, (c) Compensated Nusselt number versus Rayleigh number, (d) Compensated Reynolds 
number versus Rayleigh number. The power laws given in panels (a-b) have been derived from 
a best fit to the cases of Table 1 with Ra > 10®. In panels (c-d), the dashed lines correspond to 
the numerical solution of (6.1) and the solid lines to the numerical solution of (6.3). 


contributions by the global scalings provides a much better fit to eu and ex (figure 22). 
We can thus construct another set of equations based on the scaling laws for the total 
dissipation rates (5.7): 


-=0.248 + 7.084 

- ——^Nu =0.004 i?e +0.453 

1 + 77 + 77^ 


(6.3) 


This system of equation is once again numerically integrated to derive the scaling laws 
for Nu{Ra) and Re{Ra). The solid black lines displayed in figure 24(c-(i) and figure 25 
show that these scaling laws fall now much closer to the data. They accurately reproduce 
both Nu{Ra) and Re{Ra) for the whole range of Rayleigh numbers and the gradual 
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Figure 25. (a) versus Ra. (b) PeS versus Ra. The dashed black line correspond to the 
solution of (6.1), while the solid black line correspond to the solution of (6.3). 

change in the slopes cteg and /3eff are also correctly captured. The improvement of the 
fit to the data when using (6.3) instead of (6.1) is the direct consequence of the better 
description of the total dissipation rates by the global scalings (5.7) rather than by the 
sum of the individual scalings (5.8). 

Figures 2A{c-d) and 25 also show an extrapolation of the scaling laws, solution of (6.3), 
up to Ra = 10^^. Interestingly, the scaling laws predict that Oeff would become steeper 
than 1/3 for Ra > 5 x 10®. This transition point to an enhanced heat transport efficiency 
would then occur at much lower Ra than in RB convection in cartesian or cylindrical 
cells. For instance, the experiments by Roche et al. (2010) showed an enhanced scaling of 
Nu ^ Ra^'^^ for Ra > 7 x 10^^. The predicted effective exponent aes however seems to 
increase slightly faster than suggested by our numerical data. The possible transition to 
an enhanced heat transport regime at lower Ra than in planar geometry thus remains an 
open issue. Furthermore, the extrapolation of the obtained scaling laws to high Rayleigh 
numbers is debatable since the underlying decomposition (5.7) relies on the assumption 
of laminar boundary layers, which will not hold beyond the transition point. Future RB 
models in spherical shells that will possibly reach Ra ~ 10^® could certainly help to 
confirm this trend and check the robustness of the best-fit coefficients obtained in (5.7) 
(Stevens et al. 2013). 

7. Conclusion and outlooks 

We have studied Rayleigh-Benard (RB) convection in spherical shells for Rayleigh 
numbers up to 10® and Prandtl number unity. Because of both curvature and radial 
variations of buoyancy, convection in spherical shells exhibits asymmetric boundary lay¬ 
ers. To better characterise this asymmetry, we have conducted a systematic parameter 
study, varying both the radius ratio and the radial distribution of gravity. Two theories 
were developed in the past to determine this boundary layer asymmetry. The first one 
by Jarvis (1993) and Vangelov & Jarvis (1994) hypothesises that both boundary layers 
adjust their thickness to maintain the same critical boundary layer Rayleigh number; 
while the second one by Wu & Libchaber (1991) assumes that the thermal fluctuations 
at mid-depth are statistically symmetrically distributed. Both theories however yield 
scaling laws in poor agreement with our numerical simulations. On the contrary, we 
found that the average plume density, or equivalently the average inter-plume spacing, 
is comparable for both boundary layers. An estimation of the average plume density 
at both spherical bounding surfaces has allowed us to accurately predict the boundary 
layer asymmetry and the mean bulk temperature for the wide range of spherical shell 
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configurations explored here (77 = r^/ro spanning the range 0.2 < rj < 0.95 and gravity 
profiles g S [r/r^, 1 , (ro/r)^ {ro/rf])- 

Because of the lack of experiments and numerical models of non-rotating convection 
in spherical shells at finite Prandtl numbers, the scaling properties of the Nusselt and 
the Reynolds numbers are poorly characterised in this geometry. To further address this 
question, we have conducted numerical models in spherical shells with ry = 0.6 up to 
Ra = 10®. We have adopted a gravity profile of the form g = (co/r)^, which has allowed 
us to conduct a full dissipation analysis. One of the aims of this study was to check 
the applicability of the scaling theory by Grossmann & Lohse (2000, 2004) (GL) to con¬ 
vection in spherical shells. One of the prerequisites of this theory is the assumption of 
Prandtl-Blasius-type (PB) boundary layers. We have thus studied the temperature and 
horizontal velocity boundary layer profiles. In agreement with the previous findings by 
Zhou & Xia (2010a), the boundary layer profiles have been found to be in fair agreement 
with the PB profiles, provided the numerical simulations are analysed in a dynamical 
frame that incorporates the time and spatial variations of the boundary layers. Following 
the GL central idea, we have then decomposed the viscous and thermal dissipation rates 
into contributions coming from the fluid bulk and from the boundary layer regions. The 
detailed analysis of the individual contributions to the viscous and thermal dissipation 
rates reveals some noticeable discrepancies to the GL theory (e^ ~ ^ 

^ and ~ i?e®'®^). The total dissipation rates, however, can nevertheless be 

nicely described by the sum of bulk and boundary layer contributions that follow the 
predicted GL exponents {eu ~ aRe^ -I- bR^!"^ and tr aRe^ 6 i?e^/®). This strongly 
suggests that the inaccurate separation of the boundary layer and bulk dynamics is the 
reason for the inferior fitting of the individual contributions. These scaling laws have 
finally been employed to study the scaling properties of the Nusselt and the Reynolds 
numbers and provide laws that accurately fit the data. Although these laws exhibit a sim¬ 
ilar behaviour than experiments and numerical simulations of RB convection in cartesian 
or cylindrical coordinates; some distinction to classical RB cells have also been reported. 
Our scaling laws predict a continuous increase of the local effective slope of Nu{Ra) 
from 0.28 at Ra — 10® to 0.32 at Ra — 10® and suggest a possible enhanced heat transfer 
scaling with an effective exponent steeper than 1/3 for Ra > 5 x 10®. Similar transi¬ 
tions have been observed in some experiments, though at significantly higher Rayleigh 
numbers {Ra ~ 10^^ — 10^®, see Roche et al. 2010). 

To explore whether the spherical shell geometry is responsible for this difference, addi¬ 
tional numerical simulations at higher Rayleigh numbers are required. Ongoing improve¬ 
ments of pseudo-spectral codes for modelling convection in three dimensional spherical 
shells might help to reach spatial resolutions of the order (A/, x imax = 2048 x 2048) 
in the coming years (e.g. Schaeffer 2013). Assuming that the minimum admissible mesh 
size h has to be smaller than the global Kolmogorov scale (Grotzbach 1983; Shishkina 
et al. 2010 ) yields 

_ ^' i/i / X -(- jy -(-7^2 \ 1/4 P7.1/2 

<m - ^ J i?ai/4 (Nu-1)G4 ■ 

An extrapolation of the spatial resolutions employed in this study (Table 1) then implies 
that typical resolutions (Nr x £max = 2048 x 2048) might be sufficient to reach Ra ~ 10^® 
for the configuration we considered here (g = 0.6, g = (co/r)®). This additional decade 
in Ra might already be sufficient to ascertain the derived asymptotic scalings. 

We thank Andreas Tilgner for fruitful discussions. All the computations have been 
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carried out on the GWDG computer facilities in Gottingen and on the IBM iDataPlex 
HPG System Hydra at the MPG Rechenzentrum Garching. TG is supported by the 
Special Priority Program 1488 (PlanetMag, www.planetmag.de) of the German Science 
Foundation. 


Appendix A. Table of results for the numerical models with different 
geometry and gravity profiles 


11 


Ra 

Nu 

Re 


^h/^u 


^ui%) 

Nr ^ ^max 







g = r 

Iro 




0.2 

1 

X 

10 ® 

8.27 

723.1 

0.024/0.023 

0.015/0.021 

0.27 

0.73 

97 X 170 

0.2 

3 

X 

10 ® 

11.23 

1252.4 

0.017/0.018 

0.012/0.017 

0.29 

0.76 

129 X 341 

0.25 

2 

X 

10 " 

6.92 

407.8 

0.035/0.034 

0.022/0.026 

0.27 

0.72 

65 X 128 

0.3 

3 

X 

10 ® 

5.18 

188.1 

0.054/0.054 

0.030/0.035 

0.24 

0.69 

65 X 128 

0.3 

5 

X 

10 ® 

5.87 

242.1 

0.048/0.048 

0.027/0.031 

0.26 

0.71 

65 X 128 

0.3 

7 

X 

10 ® 

6.40 

287.2 

0.044/0.045 

0.025/0.030 

0.26 

0.72 

73 X 133 

0.3 

3 

X 

10 " 

9.38 

595.5 

0.029/0.032 

0.019/0.023 

0.28 

0.76 

73 X 133 

0.3 

3 

X 

10 ® 

18.08 

1824.8 

0.015/0.017 

0.011/0.015 

0.29 

0.81 

97 X 341 

0.35 

5 

X 

10 ® 

6.74 

274.1 

0.047/0.048 

0.027/0.031 

0.26 

0.74 

65 X 128 

0.35 

3 

X 

10 ® 

21.23 

2016.0 

0.015/0.016 

0.011/0.014 

0.30 

0.82 

129 X 341 

0.4 

1 

X 

10 ® 

5.17 

139.7 

0.067/0.071 

0.036/0.041 

0.23 

0.68 

65 X 128 

0.4 

3 

X 

10 ® 

6.72 

235.3 

0.052/0.055 

0.029/0.033 

0.23 

0.75 

65 X 128 

0.4 

5 

X 

10 ® 

7.65 

302.4 

0.045/0.048 

0.026/0.030 

0.26 

0.76 

73 X 133 

0.45 

2 

X 

10 ® 

6.71 

215.6 

0.056/0.059 

0.032/0.036 

0.24 

0.73 

65 X 128 

0.5 

1 

X 

10 ® 

6.01 

162.4 

0.066/0.070 

0.037/0.040 

0.23 

0.72 

65 X 128 

0.5 

2 

X 

10 ® 

7.22 

229.4 

0.055/0.058 

0.032/0.036 

0.24 

0.75 

65 X 128 

0.5 

5 

X 

10 ® 

9.24 

359.9 

0.043/0.045 

0.027/0.029 

0.25 

0.78 

81 X 133 

0.55 

2 

X 

10 ® 

7.71 

240.2 

0.054/0.057 

0.031/0.033 

0.24 

0.76 

65 X 128 

0.6 

1 

X 

10 ® 

6.80 

179.2 

0.065/0.067 

0.036/0.038 

0.22 

0.74 

65 X 128 

0.6 

5 

X 

10 ® 

10.55 

392.3 

0.042/0.043 

0.026/0.027 

0.24 

0.80 

81 X 133 

0.6 

5 

X 

10 ® 

10.56 

392.2 

0.042/0.043 

0.026/0.027 

0.24 

0.80 

81 X 133 

0.65 

1 

X 

10 ® 

7.11 

186.6 

0.064/0.066 

0.035/0.037 

0.23 

0.76 

73 X 170 

0.7 

7 

X 

10 ® 

6.72 

162.0 

0.070/0.072 

0.037/0.039 

0.21 

0.74 

73 X 170 

0.7 

1 

X 

10 ® 

7.41 

193.1 

0.063/0.065 

0.035/0.036 

0.23 

0.76 

73 X 170 

0.75 

1 

X 

10 ® 

7.64 

198.9 

0.063/0.064 

0.034/0.036 

0.22 

0.76 

97 X 213 

0.8 

3 

X 

10 ® 

10.60 

347.8 

0.046/0.047 

0.028/0.028 

0.23 

0.80 

97 X 426 

0.8 

4 

X 

10 " 

22.16 

1198.3 

0 .022/0.022 

0.016/0.016 

0.26 

0.85 

129 X 1024 

0.85 

7 

X 

10 ® 

7.26 

175.5 

0.068/0.069 

0.037/0.038 

0.20 

0.76 

97 X 341 

0.9 

5 

X 

10 ® 

6.73 

151.4 

0.074/0.075 

0.040/0.040 

0.22 

0.74 

97 X 426 

0.2 

1 

X 

10 ® 

11.95 

1082.7 

9 = 

0.016/0.025 

1 

0 .011/0.021 

0.29 

0.78 

97 X 170 

0.25 

2 

X 

10 " 

9.35 

572.0 

0.025/0.036 

0.016/0.027 

0.25 

0.78 

81 X 133 

0.3 

7 

X 

10 ® 

8.15 

377.8 

0.033/0.046 

0.020/0.031 

0.26 

0.78 

73 X 133 

0.35 

6 

X 

10 ® 

8.79 

383.9 

0.035/0.046 

0.021/0.030 

0.26 

0.79 

73 X 133 

0.4 

3 

X 

10 ® 

7.98 

288.7 

0.042/0.054 

0.025/0.032 

0.24 

0.79 

65 X 128 

0.4 

5 

X 

10 ® 

9.19 

371.4 

0.036/0.047 

0.022/0.029 

0.27 

0.80 

65 X 128 
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0.45 

5 X 10® 

9.95 

395.4 

0.036/0.046 

0.022/0.029 

0.27 

0.80 

73 X 128 

0.5 

5 X 10® 

10.57 

421.7 

0.036/0.045 

0.023/0.029 

0.26 

0.80 

81 X 133 

0.55 

5 X 10® 

11.09 

434.8 

0.037/0.044 

0.023/0.028 

0.26 

0.81 

81 X 133 

0.6 

3 X 10® 

10.04 

343.3 

0.042/0.049 

0.026/0.030 

0.24 

0.80 

81 X 133 

0.6 

5 X 10® 

11.68 

442.7 

0.036/0.042 

0.023/0.027 

0.26 

0.81 

81 X 133 

0.65 

3 X 10® 

10.48 

355.3 

0.042/0.048 

0.026/0.029 

0.24 

0.80 

81 X 133 

0.7 

1 X 10® 

7.85 

210.0 

0.058/0.064 

0.033/0.036 

0.22 

0.77 

73 X 213 

0.75 

3 X 10® 

10.82 

363.0 

0.043/0.047 

0.026/0.029 

0.24 

0.80 

97 X 341 

0.8 

3 X 10® 

10.98 

367.4 

0.043/0.046 

0.027/0.028 

0.24 

0.80 

97 X 426 

0.85 

1 X 10® 

8.22 

217.5 

0.059/0.062 

0.033/0.035 

0.21 

0.76 

97 X 426 

0.9 

5 X 10® 

6.83 

155.2 

0.072/0.074 

0.039/0.040 

0.21 

0.74 

97 X 426 





9 = {i"c 





0.2 

5 X 10® 

5.85 

171.9 

0.031/0.089 

0.017/0.047 

0.25 

0.80 

49 X 85 

0.2 

7 X 10® 

6.45 

206.2 

0.028/0.082 

0.016/0.046 

0.24 

0.80 

49 X 85 

0.2 

1 X 10® 

7.17 

248.8 

0.026/0.076 

0.014/0.044 

0.24 

0.80 

61 X 106 

0.2 

1.5 X 10® 

8.03 

308.1 

0.022/0.068 

0.012/0.041 

0.28 

0.81 

81 X 170 

0.2 

2 X 10® 

8.78 

351.9 

0.020/0.064 

0.011/0.039 

0.26 

0.82 

81 X 170 

0.2 

1 X 10^ 

14.19 

794.2 

0.013/0.040 

0.008/0.027 

0.30 

0.84 

97 X 256 

0.2 

3 X 10^ 

20.25 

1369.5 

0.009/0.029 

0.006/0.021 

0.31 

0.84 

97 X 341 

0.25 

5 X 10® 

6.29 

180.2 

0.035/0.087 

0.019/0.047 

0.26 

0.79 

49 X 85 

0.3 

5 X 10® 

6.60 

185.7 

0.038/0.085 

0.021/0.046 

0.22 

0.80 

65 X 128 

0.3 

1 X 10® 

8.08 

263.2 

0.031/0.071 

0.017/0.041 

0.23 

0.81 

65 X 128 

0.3 

3 X 10® 

11.11 

458.1 

0.022/0.053 

0.013/0.034 

0.29 

0.83 

73 X 133 

0.3 

1 X 10^ 

16.07 

803.7 

0.015/0.037 

0.010/0.025 

0.28 

0.86 

97 X 256 

0.3 

5 X 10^ 

25.68 

1810.7 

0.010/0.024 

0.007/0.018 

0.30 

0.86 

161 X 512 

0.35 

5 X 10® 

6.87 

187.4 

0.041/0.083 

0.023/0.045 

0.23 

0.79 

65 X 128 

0.4 

5 X 10® 

7.15 

189.0 

0.044/0.081 

0.024/0.043 

0.24 

0.78 

65 X 128 

0.4 

8 X 10® 

8.07 

237.7 

0.038/0.072 

0.022/0.041 

0.22 

0.80 

65 X 128 

0.4 

1 X 10® 

8.59 

266.6 

0.036/0.068 

0.021/0.039 

0.26 

0.81 

65 X 128 

0.4 

3 X 10® 

11.72 

459.1 

0.026/0.051 

0.016/0.032 

0.27 

0.83 

97 X 170 

0.45 

5 X 10® 

7.26 

192.7 

0.046/0.080 

0.026/0.044 

0.24 

0.77 

65 X 128 

0.45 

7 X 10® 

7.98 

226.9 

0.042/0.073 

0.024/0.041 

0.23 

0.78 

65 X 128 

0.5 

5 X 10® 

7.29 

190.6 

0.049/0.079 

0.028/0.044 

0.23 

0.76 

61 X 106 

0.5 

7 X 10® 

8.01 

225.7 

0.045/0.072 

0.026/0.041 

0.23 

0.78 

65 X 128 

0.5 

1 X 10® 

8.83 

270.1 

0.041/0.066 

0.024/0.039 

0.23 

0.79 

81 X 170 

0.55 

5 X 10® 

7.31 

187.3 

0.052/0.078 

0.029/0.043 

0.24 

0.78 

61 X 106 

0.65 

5 X 10® 

7.36 

182.2 

0.057/0.076 

0.031/0.042 

0.21 

0.77 

61 X 106 

0.7 

1.5 X 10® 

5.22 

97.7 

0.083/0.106 

0.042/0.053 

0.19 

0.72 

73 X 170 

0.7 

3 X 10® 

6.35 

138.2 

0.068/0.087 

0.036/0.046 

0.20 

0.75 

73 X 170 

0.7 

5 X 10® 

7.32 

178.3 

0.059/0.075 

0.032/0.041 

0.21 

0.76 

81 X 266 

0.7 

1 X 10® 

8.83 

251.0 

0.049/0.063 

0.029/0.036 

0.22 

0.79 

73 X 213 

0.75 

5 X 10® 

7.26 

174.9 

0.062/0.075 

0.033/0.041 

0.22 

0.75 

81 X 266 

0.8 

5 X 10® 

7.20 

171.2 

0.064/0.074 

0.035/0.040 

0.20 

0.75 

81 X 266 

0.8 

1 X 10® 

8.71 

240.8 

0.053/0.062 

0.030/0.035 

0.22 

0.78 

97 X 341 

0.8 

3 X 10® 

11.81 

409.2 

0.039/0.046 

0.024/0.028 

0.23 

0.81 

97 X 426 

0.85 

1 X 10® 

8.62 

235.5 

0.055/0.061 

0.031/0.034 

0.22 

0.78 

97 X 426 

0.9 

1 X 10® 

4.47 

72.2 

0.108/0.116 

0.053/0.056 

0.17 

0.68 

97 X 266 

0.9 

5 X 10® 

7.04 

163.8 

0.069/0.074 

0.037/0.040 

0.21 

0.75 

97 X 426 

0.9 

1 X 10® 

8.54 

230.2 

0.057/0.061 

0.032/0.034 

0.22 

0.78 

97 X 512 
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0.95 
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1 

X 
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5 

X 

10 ’^ 

13.04 

0.45 

3 
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10.25 

0.5 

1 

X 

10 ’^ 

6.95 

0.5 

3 

X 
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9.53 

0.55 

3 

X 

10 ^ 

8.88 

0.6 

3 

X 

10 ^ 

8.41 

0.6 

5 

X 

10 ^ 

9.69 

0.6 

7 

X 

10 ^ 

10.69 

0.6 

1 

X 

10 " 

23.28 

0.65 

3 

X 

10 ’^ 

7.97 

0.7 

3 

X 

10 ’^ 

7.60 

0.7 

5 

X 

10 ’^ 

8.76 

0.7 

7 

X 

10 ’^ 

9.66 

0.75 

5 

X 

10 ^ 

8.34 

0.8 

7 

X 

10 ^ 

8.80 

0.8 

7 

X 

10 ^ 

8.78 

0.85 

7 

X 

10 ^ 

8.41 

0.9 

1 

X 

10 ® 

8.95 


70.4 0.112/0.116 

9 = {r, 

206.2 0.017/0.117 

265.7 0.015/0.100 

313.9 0.013/0.089 

238.3 0.020/0.105 

212.7 0.025/0.107 

251.9 0.023/0.095 

440.1 0.015/0.070 

1327.6 0.008/0.035 

307.3 0.023/0.083 

189.5 0.034/0.104 

327.1 0.024/0.076 
422.0 0.020/0.067 
290.0 0.029/0.077 

153.4 0.046/0.106 

265.8 0.033/0.079 

239.3 0.038/0.081 

216.6 0.043/0.081 

278.4 0.037/0.071 

329.1 0.033/0.065 

1199.9 0.015/0.031 

199.8 0.047/0.082 

183.9 0.052/0.082 

236.1 0.045/0.072 

279.4 0.041/0.066 

218.8 0.050/0.072 
240.0 0.049/0.066 

239.7 0.049/0.066 

223.2 0.054/0.066 

248.9 0.052/0.060 


0.054/0.056 

0.17 

./rf 


0.009/0.055 

0.24 

0.007/0.048 

0.23 

0.007/0.043 

0.28 

0.011/0.055 

0.26 

0.013/0.056 

0.24 

0.012/0.050 

0.23 

0.009/0.040 

0.24 

0.005/0.023 

0.29 

0.013/0.045 

0.23 

0.018/0.055 

0.23 

0.014/0.043 

0.23 

0.012/0.039 

0.27 

0.016/0.042 

0.26 

0.025/0.053 

0.22 

0.019/0.044 

0.23 

0.022/0.045 

0.25 

0.024/0.044 

0.23 

0.021/0.040 

0.24 

0.020/0.038 

0.24 

0 .011/0.022 

0.28 

0.026/0.044 

0.24 

0.029/0.044 

0.22 

0.026/0.040 

0.22 

0.024/0.038 

0.22 

0.028/0.040 

0.21 

0.028/0.038 

0.22 

0.028/0.037 

0.23 

0.030/0.037 

0.22 

0.030/0.034 

0.21 


0.68 

97 

X 

512 

0.85 

65 

X 

128 

0.85 

65 

X 

128 

0.85 

65 

X 

133 

0.82 

65 

X 

128 

0.83 

65 

X 

128 

0.83 

65 

X 

133 

0.84 

65 

X 

133 

0.86 

161 

X 

: 426 

0.83 

65 

X 

133 

0.80 

65 

X 

133 

0.83 

65 

X 

133 

0.83 

65 

X 

133 

0.81 

65 

X 

133 

0.76 

65 

X 

133 

0.80 

65 

X 

133 

0.78 

65 

X 

133 

0.78 

65 

X 

133 

0.80 

73 

X 

170 

0.81 

73 

X 

170 

0.84 

129 

1 X 

: 341 

0.77 

65 

X 

170 

0.77 

65 

X 

170 

0.78 

65 

X 

170 

0.78 

65 

X 

170 

0.78 

65 

X 

256 

0.79 

65 

X 

256 

0.78 

97 

X 

426 

0.78 

97 

X 

426 

0.78 

97 

X 

512 


Table 2: Summary table of Pr = 1 numerical simulations with 
various radius ratio ry and gravity profiles g{r). 
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